NeuralFluid: Nueral Fluidic System Design and Control with Differentiable Simulation
Authors: Yifei Li
Yuchen Sun
Pingchuan Ma
Eftychios Sifakis
Tao Du
Bo Zhu
Wojciech Matusik
Venue: NeurIPS 2024 poster (2024)
Keywords: differentiable simulation
fluid learning
fluid simulation
fluid control
fluid optimization
fluid design
inverse problems
navier-stokes
adjoint method
fluid dynamics
Primary area: machine_learning_for_physical_sciences

Abstract: We present NeuralFluid, a novel framework to explore neural control and design of complex fluidic systems with dynamic solid boundaries. Our system features a fast differentiable Navier-Stokes solver with solid-fluid interface handling, a low-dimensional differentiable parametric geometry representation, a control-shape co-design algorithm, and gym-like simulation environments to facilitate various fluidic control design applications. Additionally, we present a benchmark of design, control, and learning tasks on high-fidelity, high-resolution dynamic fluid environments that pose challenges for existing differentiable fluid simulators. These tasks include designing the control of artificial hearts, identifying robotic end-effector shapes, and controlling a fluid gate. By seamlessly incorporating our differentiable fluid simulator into a learning framework, we demonstrate successful design, control, and learning results that surpass gradient-free solutions in these benchmark tasks.

NeuralFluid: : Nueral Fluidic System Design and Control with Differentiable Simulation
Yifei Li MIT CSAIL
Yuchen Sun
Pingchuan Ma MIT CSAIL
Georgia Institute of Technology
Eftychios Sifakis
Tao Du
University Wisconsin-Madian
Tsinghua University, Shanghai Qi Zhi Institute
Bo Zhu Georgia Institute of Technology
Wojeiech Matusik MIT CSAIL
Abstract
We present NeuralFluid novel framework to explore neural control and design complex fluidic systems with dynamic solid boundaries. Our system features a fast differentiable Navier-Stokes solver with solid-fluid interface handling, a low- dimensional differentiable parametric geometry representation, a control-shape co-design algorithm, and gym-like simulation environments to facilitate various fluidic control design applications. Additionally, we present a benchmark design, control, and learning tasks high-fidelity, high-resolution dynamic fluid environments that pose challenges for existing differentiable fluid simulators. These tasks include designing the control of artificial hearts, identifying robotic end- effector shapes, and controlling a fluid gate. By seamlessly incorporating our differentiable fluid simulator into learning framework, we demonstrate successful design, control, and learning results that surpass gradient-free solutions in these benchmark tasks.
1 Introduction
Complex fluidic systems play an important role many engineering and scientific disciplines, en- compassing applications at different length scales ranging from biomedical implants [1, microfluidic devices (2), hydraulic devices to and flying robots 3. Understanding these fluid-solid coupling mechanisms in nature and mimicking their control strategies in artificial designs is essential for advancing our control and design capabilities to synthesize novel solid-fluid systems.
Devising neural control algorithms to accurately manipulate the behavior of complex fluidic system and optimize performance remains challenging due to the intricate interplay between device geometry, control policies, flow dynamics, and the inherent physical and optimization constraints unique each fluidic system. On one hand, differentiable simulation fluid-systen interactions are inherently difficult because simulation is dynamic, involving sequence of forward and backward steps interleaved with control signals that are computationally expensive. the other hand, naively employing traditional control algorithms, mainly derived from their solid counterparts, to control fluidic systems remains difficult due to characterizing the infinite degrees freedom fluid flows and their interactions with solid boundaries. The co-design of fluid-solid systems, involving both shape and control, is critical to exploring the optimal performance of these systems.
Currently, the machine learning community lacks computational Gym-like 4 environment to facilitate the exploration of fluidic systems manifesting strong solid-fluid interactions and controllable
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

Geometry & Neural Control
Differentiable Navier-Stokes Simulation
Loss Computation
Simulated
M
Target
Controller
Geometry
Optiised
Forward Simulation
Backward Propagation
Time
Figure Pipeline Overview. (1) Our pipeline starts with an initial parametric geometry and a neural network parameterized controller. (2) The fluid dynamics is then simulated using a dynamic Navier-Stokes solver. (3) The performance of the design and control is evaluated using loss function, the gradients of which are then prok-proaggated through
end-to-end differentiable framework. (4) The gradient- based optimization iteratively improves the geometry and control to achieve the task goal. This pipeline allows for efficient geometry and control co-optimization.
dynamic boundaries. Recent literature in robotic learning (e.g.. (5) has established unified multi- physics differentiable simulation platforms to facilitate learning control policies for various fluid interactions in daily scenarios. Similar ideas can be observed in 6 7 8) where differentiable simulation plays a central role in accommodating various design and optimization tasks of dynamic systems involving fluid dynamics. However, despite these inspiring advances, learning the control policies and exploring the optimal performance of a dynamic fluidic system with complex boundary conditions remains difficult due to their inherent complexities in differentiating solid boundary behaviors and optimizing their fluidic consequences due to these boundary motions.
This paper presents a novel framework for fully automated pipeline aimed at devising neural controls for complex fluidic systems with dynamic boundaries. Our framework is designed to robustly control complex fluidic systems that consist of externally driven soft boundaries and internal complex flow behaviors, such as those systems underpinning an artificial heart or a microfluidic device.
NeuralFluid consists of three critical components to enable neural control of a complex fluidic system. First, we devise a differentiable geometry representation to offer an expressive design space while remaining low-dimensional, enabling efficient exploration by the optimization algorithm. Second, we implement a differentiable fluid simulator with solid-fluid interface handling to accurately characterize the dynamic fluid behavior and predict its spatiotemporal impact on the moving boundaries. We back-propagate gradients at the solid-fluid interface to extend gradient computation to the geometry iso-surface. Last, we provided an optimization framework to efficiently search the design space. considering the underlying fluid dynamics and boundary conditions.
Our pipeline features niw-dimensional parametric geometry representation capable of expressing complex shapes and a differentiable Navier-Stokes simulator with geometry gradient computation for predicting dynamic fluid behavior in response control signals. In addition, our pipeline leverages gradient-based optimization for efficient design space exploration, co-optimization of the device geometry and control, and accurate performance evaluation of the design under dynamic flows. To showease the practical implications and versatility of our approach, we have established a suite of Gym-like [4] environments. These benchmarks are designed to test applications in robotics and engineering, facilitating advancements in system identification, optimization of end-effector shapes and controls, and the dynamic optimization of structures such as artificial hearts within a closed-loop control framework. We showcase the effectiveness of our pipeline in facilitating different design and control tasks, including amplifier, fluidic switch, flow modulator, shape and position identification, closed-loop control of water gate and artificial heart.
We summarize main contributions follows:
Development of a fast differentiable Navior-Stokes simulator for optimization in 2D and 3D scenes.

---

Development of low- -ww-dmenssionl differentiable parametric geometry representation for complex shapes embedded into the differentiable simulation pipeline.
Gradient computation extension to geometry iso-surface enable control and geometry co-design and iso-surface optimization.
Gym-like [4] environments and benchmarks to demonstrate applications in robotics and
engineering, including the design of amplifier, fluidic switch, flow modulator, geometry system identification, and closed-loop control of fluid gate and artificial heart.
2 Method
2.1 Pipeline Overview
We present overview of our method in Fig. D Our pipeline defines the designs with low- dimensional parametric geometry representation (Sec. (22) The behavior and performance of the design in the fluid environment (Sec. 2.3) Both components are embedded in a gradient- -hased optimization framework that co- optimizes both the geometric design and the control signal until convergence.
evaluated by dynamic differentiable Navier-Stokes simulator
2.2 Geometry Representation
We represent our geometry with a low-dimension rep- resentation. Take the illus- tration the inset figure as an example, here we intro- duce the representation on high level, and refer the readers to the appendix for the full details. We param- eterize closed 2D surface using its center and set of connected Bezier curves with their control points de- fine in polar coordinates P; for E curve spanning 25 radians in the polar coordinate system. This representation offers compact way of defining diverse geometries. We further parameterize closed 3D surface using a list of 2D surfaces defining the key cross-sections of the geometry along extrusion axis of the local object frame, where each 2D surface is parameterized as described above. The parameterization includes
Co
P1
Pi
(c,, c,)
a) Closed 2D Surface
b) Closed 3D Surface
3D Geometry
where every two control points define a 2D Bezier
zo and 21. which determines the Z plane of the first and last cross-section, along with the parameters for each key 2D cross-section, which are assumed to be evenly spaced between E [30, 21] The continuous geometry interpolates the key cross-sections along the z-axis. Finally, we construct more complex 3D geometries using operations from Constructive Solid Geometry (CSG): Union and intersection, which allows define a 3D parametric heart model using the union of four sub-geometries.
2.3 Differentiable Navier-Stokes Simulation
Our fluid dynamics is governed by the incompressible Navier-Stokes equations. These consist the momentum equation (Eq. la), accounting for temporal changes in velocity (u), advective acceleration, viscous dissipation, and pressure (p) gradient forces for an incompressible fluid with fluid density and kinematic viscosity The incompressibility condition (Eq. [bb requires the divergence of the velocity field must be zero to enforce the conservation of mass:
du
(u V)u v Vp.
(1a)
Øt
v.u=0
(1b)
3

---

2.3.1 Numerical Simulation
We build the fluid simulator by leveraging the operator-spliting method [0] A single simulation step comprises three sub-steps: advection, viscosity, and projection. See the Appendix for details on time discretization. The simulation domain is discretized grid [1I], with pressures stored at cell centers and velocities at cell faces. By employing the finite- difference scheme the MAC grid cells and faces, construct the matrix Ae G for gradient operator and its negative transpose - Az G for divergence operator. In the following sections, capital letters will refer to matrices or the flattened vectors induced by the fields denoted by the corresponding lowercase letters in Appendix B.
standard Marker and--ad-aa Ced- (MAC)
Advection We employ the semi-Lagrangian advection scheme, where the advected velocity field U Ü is linear interpolation of the velocity field U". The interpolation position function is function of U can be put into a matrix form B, which results in:
n+1 B(U")
Viscosity For incompressible fluid with constant viscosity coefficient, the viscous force density is the product of the Laplacian of velocity and the viscosity coefficient. For each axis, the Laplacian of the corresponding velocity component is calculated on grid points using the finite difference method:
[
v v G G
ūnt! ))
(3)
Ax
Projection The projection step ensures the incompressibility of the fluid. Solid un- aligned with the grid may intersect with grid faces, which can be captured with a cut-cell method [12]. We introduce "+1 to represent the fluid proportion
B
1.80
grid face. The solid's signed distance function (SDF) dn+1 and the velocity " ++ can derived from the solid geom- etry. We first use marching cube compute the geometry zero contour. Next, we identify the intersection points between this contour and the grid
2 j.k)
D2
faces and compute a" n+1 based Φ. For instance, in the inset figure, grid face (i 2. k)) is cut by
the contour, then 14
SABE SABCD
[AE]
LAF
2 IAB
IAD
o
7+1 The volume change rates for fluid and solid at grid cell (i,j, k). denoted as Yfask and Teijk respectively, equal to the sum of flux on the cell's surrounding faces, which can be calculated using an+1, fluid velocity u=+1, and solid velocity u
n+1 The incompressibiliity condition gives requires the sum of Tfijj and n+1 Ta
to be zero, which gives
Al AAr
- Gpa+1 G n1
+ G Snt S
where sh+1 is diagonal matrix induced by a n+1 and P is the pressure. After solving the linear system, the fluid velocity is updated based on the pressure values:
Al
Gp
(5)
pAr
2.3.2 Back-propagatiion through Time
We construct our back buck-propagation algorithm mirror the sequence of operations carried out in the forward pass but in reversed order.
Given the gradients of the loss function J with respect to the velocity field u at time step denoted by - our goal is to compute the corresponding gradients at time step ~U
2J

---

Projection
We begin by reversing the projection step to back-pue proppapaeaaeee
aJ at
to derive
aj ape and
a Back-propagating through Eq
2j
gives
OJ
At
DJ
G. páv a0
(6)
0p
w.r.t U n+1
We can pack-propagate the adjoint of Eq.
by defining the adjoint variable y and derive
0J
0J
yyG
aÜ
n+1
au
where A AL L Sn Sg = G S 1
n+
+
and y is computed by solving
pår
the linear system Ayj (pp-T)
Viscosity and Advection Back propagating through viscosity and advection simply involves back- through Eq. and Eq. which allows us to derive:
propagating at"+
8J
0J
OJ
1
vit Ar
ŌB T
G²
U + B(U))).
aU
aU 3
The above equations provide the outline of the back-propagation process through single time step of from time step n + n. To compute the gradients of the loss function any time step, we iterate the back precppropagttion process over the full sequence of time steps.
2.3.3 Back-propagation through Geometry
The parametric geometry affects simulation through the solid-fluid boundary during the projection step in Eq. 4 Specifically, the SDF of the geometry d=+1 affects the volume matrix s and the velocity (in the case of moving geometry) of the geometry U affects the boundary condition. We can back-propagate
w.r.t these two parameters derive apn
8J
u yG (Ü* a+t
OJ
Un+1) + y
aA p=+1 as
(9a)
JJ
au yG
(9b)
87 Further back-propagating u5 first through the SDF & then through the distance computation and 00 through geometry velocity function allows us to optimize through the geometry iso-surface.
2.3.4 Neural Fluid Control
We can train
parameterized closed-loop fluid controllers with gradients fully computed at both geometry and velocity throughout time. We parameterize our controllers with a two-layer MLP. The controller takes input the observation of the fluid velocity field at each frame and outputs dynamic control signals that affect the geometry through our parametric geometry presentation, which further affects the flow field. Our fully differentiable framework allows gradient-bant based methods to train the controller efficiently. We implemented the backbone of our code in C++ and CUDA for computational efficiency. We derived gradients for the geometry and simulation module analytically, then exposed the differentiable simulation framework through pybindl integration with deep learning libraries, which in case is PyTorch [4
13 to enable seamless
3 Benchmarks and Applications
In this section, we introduce our fluidic design benchmarks and environments. We warp our environ- ments using the standard protocol the gym to facilitate learning practices. We present six fluidic design and control tasks (Fig. 2 to assess the effectiveness of our computational pipeline for fluidic
5

---

Table 1: Task Specifications. We summarize the simulation and optimization configuration for the design tasks shown in Sec. 3 and report the initial and optimized loss. We note that because our implementation adopts CFL condition for numerical stability during simulation, the actual steps simulated and back-propagated are higher than the numbers shown in "# Frames".
Resolution
# Frames
# Param.
BL DDsdgn
Bb Conttro
Loss (Ly) Optimized
Initial
Amplifier Switch
64
5
13,401
0,005
64 x 64
120
10
13.162
1.893
Shape Identifier
128 x 128
10
10
V
70.152
0.759
Flow Modulator
40 x 40 40
100
34
V
6.118
0.069
Neural Gate
x
4.5k
7.318
0,000
Neural Heart
48 x 48 48
180
7.1k
1.086
0,004
Amplifier
Switch
Shape Identifier Flow Modulator
Neural Gate
Neural Heart
Figure 2: Tasks Overview. In each task, the blue dashed line represents the inlet, the red dashed line indicates the outlet, the white arrows show the flow direction, and the orange shapes and arrows denote the geometry and its motion direction.
system design and learning. A comprehensive illustration of these design scenarios, including the visualization of the optimization process, is provided in Appendix Sec|C]and our supplemental video. Initial conditions are set for all optimizations using randomly sampled values. We use Adam as our optimizer. We summarize the simulation configuration and optimization configuration as well as relevant statistics Table
3.1 Task Overview
This design problem aims to amplify parallel horizontal inflow by three times from an initial velocity of 5 units to 15 units. The device boundaries are parameterized as two symmetrically placed cubic Bezier curves. The design variables are the control points and endpoints of the two curves. The loss function defined as the last frame L2 norm the difference between the target and optimized fluid velocity norm. We visualize the initial and optimized designs in Fig. 3h. We overlay the design and the corresponding velocity field (colored by the norm) for both iterations.
Amplifier
Shape Identifier This task provides example of system identification in fluid environment by identifying the shape and position of a geometry, given observations of the flow field. We randomly
Initial
Optimized
Velocity Streamline
Inlet wall
ex0 Aut
QUILE
objection
9.283
& V
Optinizd
Fow Concrolie- 3x
Foo Concrstle-54 Sy.
5.748
50
1000
Amplilier
Iterationn
Db Flow Modulator
Figure 3: (a) Visualization of Amplifier. (b) Visualization of Flow Modulator.

---

Monde1
Moucle
Marcle
Maxcled
Optimbed
Yarget
Velecity
150
1.00
Model
-105
10.0
5.04
0.84
Time
5.00
0.06
0,0 Time
Central
Veloctt
5.100
Monde
Mesde
1000
Figure 4: Artificial Heart. Left: visualization of the domain and the location of the muscles. Middle: Optimized control policy rollout visualization. Right: Optimization results visualization. The top and bottom diagrams visualize the cosine and the ECG target variants.
initialize the geometry in the domain. We define the loss function the sum of the L2 norm of the velocity field difference to the observed ground-truth flow field across time. The optimization successfully reconstructs the shape and its position with random initialization (See Exp. 4.2
Switch This task simultaneously optimizes the geometry and the constant rotational speed of a 2D switch, allowing dynamic regulation of outlet flow velocity. horizontal inflow at the left interacts with the switch, splitting into distinct streams towards the right side. The goal is let the time-dependent average velocity norm of the upper stream align with a predetermined flow profile. We define the loss function the of the L2 norm of the difference between the average velocity norm of the fluid and the target velocity at each frame. The optimized design and control successfully generate linearly increasing upper stream velocity norm profile, matching the specified target.
Flow Modulator This task optimizes the geometry and control of rotating 3D flow controller to achieve a target average outlet flow (kx inflow) at the domain's right boundary by the end of the simulation. The controller's geometry is parameterized by four 3D cross-scctions, with rotation controlled by a sinusoidal function. We define the loss as the L2 norm the difference between the average outlet and target velocity at the final frame. Fig. Bb illustrates the task: the top left shows initial parameters and task specification, bottom left shows optimization trajectories, and the right visualizes the optimized geometry and velocity streamline with two variants (k: 3, 5).
We learn closed-loop controller for 3D fluid gate moving horizontally. The controller parameterized two-layer MLP. observes the current outflow velocity through the gate and outputs the next frame motion offset to control the outflow velocity to match target.
Neural Gate Controller
3.2 Scalability Complex Fluid Fields: A Case Study of Artificial Heart
We showcase the scalability of our method through artificial heart design and control learning task (Fig. 4.. Artificial heart development is difficult due to the complex blood flow movement within the heart. This study provides first step in studying the heart's control strategies. We train a closed-loop controller that outputs the per-time-step contraction signal of the four muscles of a simplified heart model so that the outlet velocity matches a pre-defined target profile. The controller's states include temporal encoding of the current time step and the current outflow norm, and they are parameterized using a two-layer MLP. The heart's geometry is parameterized as the union of the two inlets, one outlet, and the heart chamber. In two variants of the task, one target flow profile is parameterized using cosine curve (Fig top)), and one target flow profile mimics the shape of an electrocardiogram (Fig 4 bottom.. We define the loss function as the sum of the L2 norm of the
difference between the average velocity norm the fluid and the target velocity at each frame. In both variants, the trained controller successfully outputs signals that generate blood flow that matches the target, demonstrating the effectiveness of our gradient-based optimization framework. We further visualize the rollouts of the trained controllers at Fig 4ig4lo left.
7

---

Table 2: Time Performance. Our method achieves one order of magnitude speedup across all
resolutions compared to PhiFlow in both forward simulation and backward gradient propagation.
Forward
Backward
Resolution
PhiFlow (s)
Ours (s)
Speedup
PhiFlow (s)
Ours (s)
Speedup
32 32
32
1.282
0.024
53.4x
1.546
0.095
16.3x
40 40 40
1.741
0.039
44.6x
1.983
0.121
16.4x
48 48 x 48
2.227
0.068
32.8x
2,412
0.158
15.3x
64 x 64 x 64
3.145
0.105
30.0x
4.094
0.301
13.6x
Table 3: Memory and time performance comparison with DiffTaichi
Memory (MB)
Forward Time (8)
Backward Time (s)
Resolution
DiffTaichi
Ours
DiffTaichi
Ours
DiffTaichi
Ours
32 x 32 32
685
292
0.081
0.024
0.074
0.027
40 x 40 40
1136
308
0.146
0.039
0.133
0.041
48 48
1805
322
0.228
0.068
0.183
0.064
64
5005
405
0.435
0.105
0.363
0.117
Experiments
4.1 Effects of Initialization on Optimization
This experiment studies the effect of random initialization in our fluid optimization tasks. Specifically, we conduct an experiment on the 3D heart controller task, which utilizes neural network with 7,100 parameters. For this task, we initialized the network parameters with five different random seeds, tracking convergence under each condition. In the accompanying figure, we plot extended version of the training trajectory, scaled logarithmically for better visualization, to compare the convergence of different random initializations. However, in practice, our method achieves the objective driven by the loss within tens of iterations, as evident from the steep initial descent in the optimization curve. As shown in Fig. 5 left, the optimization consistently converges across all seeds, despite variations in initial network behavior. This consistency indicates the robustness our method to random initialization even in high- high-dimensiooal optimization spaces, supporting its application to complex tasks in differentiable physics. These curves offer valuable guidelines for practitioners using method in their deployments: the gradients from approach robust to parameersssstrss and scalable to high-dimensional optimization problems.
4.2 Gradient-Based Gradient-Free Optimization
We study the effectiveness of our gradient-based method against gradient-free optimization methods in the Neural Heart task (Fig. 1 right). We choose Proximal Policy Optimization (PPO) [15] the baseline for reinforcement learning and CMA-ES [16] for evolution strategies. This environment is particular challenging due the sensitivity of the bloodflow to the control signal changes across time, which could result in large flow field change if adjacent rollouts have large changes. We initialize all methods output stochastic control signals of small noise for stable initial simulation. Our gradient-based method quickly converges to near zero after 60 epochs, while both gradient-free methods struggle in this environment. We argue that the rapid and successful convergence stems from the clear gradient provided by our method. Note that our differentiable optimization pipeline depends on the differentiability of the loss function (e.g., the imitation loss we used). This can be problematic if the objective is too complex to be characterized in a differentiable manner, in which case gradient-free methods are better alternatives. However, our framework will still excel due to its outstanding forward simulation speed, which we will elaborate next.

---

aa) Initialization
(b) Ours v.s. Gradient-free Methods
Seed
Seed
10'
15
Seed
10"
Seed
Loas
10
Seed
Loss
10'
10
Ours
10
10" CMA-ES
#lterations
#fterations
Figure 5: Ablation Studies. Left: Optimization trajectories for Neural Heart with 7100 parameters under different initialization. Iterations are visualized on log scale. Right: Log scaled loss-iteration curves of our gradient-baad based method and other gradient- optimization methods.
4.3 Time Performance Profiling and Comparison with PhiFlow
In this experiment, we demonstrate the performance efficiency of our framework through a compar- ison with PhiFlow. While PhiFlow operates with a TensorFlow-GPP backend, our framework is implemented in CUDA C++ and features a high-performance Geometric-Mulligee poultigrid-Precondiiionnee Conjugate- Gradient (MGPCG) Poisson solver [7]. To address the needs of differential operators and interpolations, which require access to neighboring cells in all directions, we divide the simulation domain into cubic blocks, each corresponding to a CUDA block. When launching CUDA kernel, simulation data for each block is first loaded into shared memory, allowing efficient computation directly in shared memory and reducing global memory accesses. Additionally, to increase memory throughput, each block's 's data is stored consecutively in global memory. Our matrix-free MGPCG solver has a faster convergence rate than PhiFlow's 's Conjugate Gradient solver and uses a hicrarchical grid data structure, with custom CUDA kernels for prolongation and restriction operations between coarse and fine grids.
We benchmark both the one-step forward simulation time and gradient pack-propagation time at different resolutions, as shown Table The experiment runs on workstation with an NVIDIA RTX A6000 GPU, where our framework consistently outperforms PhiFlow by an order of magnitude across all resolutions, benefiting both gradient-based and gradient-free optimization techniques. This performance improvement is particularly advantageous in fluid simulation applications. including roboties and video generation. Additionally, our system's gym protocol compatibility + makes straightforward for practitioners to integrate and test our library. We plan to release code and documentation upon acceptance.
4.4 Memory and Time Performance Profiling and Comparison with DiffTaichi
Here we compare our solver with DiffTaichi, a differentiable programming framework, to highlight the benefits of our approach in terms of scalability and efficiency. Our method, designed specifically for differentiable fluid simulation, uses manually derived gradients, avoiding the need store intermediate computational graph at each timestep, unlike DiffTaichi, which relics on automatic differentiation. Additionally, adjoint derivation for the projection solve step is independent of solver iterations, making our approach well-suited for advection-projees To demonstrate this, implemented Conjugate Gradient (CG) solver for the projection step DiffTaichi and compared time and memory performance across four grid resolutions a 3D optimization scenario (Table 3. Our results show that solver requires substantially less memory, with up to 12 times less memory usage than DiffTaichi at 64 x 64 x 64 resolution. This reduction in memory stems from eliminating the need to store intermediate values during each CG iteration, making our solver particularly suitable for high-resolution, long-term optimizations.
fluid simulations.
5 Related Work
Flow Control and Optimization Beginning with the pioneering work of [18, a vast literature has been devoted to the optimization of fluid systems [19]. Given predefined design domain with

---

boundary conditions, typical optimization objective is to maximize some performance functional of fluid system (e.g,, the power loss of the system) constrained by the physical equations. Similar to a conventional structural optimization problem, the design domain is discretized. The optimization algorithm decides for each element whether it should be fluid or solid to optimize some performance functions such as power loss. Examples of flow optimization applications include Stokes flow
[8 20 22). steady-state flow [23], weakly compressible flow [24], unsteady flow [25], channel
flow [26]. ducted flow (27], viscous flow (28). fluid-structure interaction (FSI) [29] 30 31, fluid-
thermal interaction 32 [d
microfluidics [34], aeronautics [35 36., and aerodynamics [37 38],
to name few. 39 developed dynamic differentiable fluid simulator and integrated the pipeline with neural networks for learning controllers. In computer graphics, [40] developed a differentiable framework to simulate and optimize flow systems governed by design specifications with different types of boundary conditions, while [8] developed an anisotropic material model to handle different boundary conditions using topology optimization framework. Both systems focus the Stokes flow model and have not explored applications with dynamic flow system. 41 adapted the adjoint method to control free-surface liquids.
Differentiable Physics Simulation Differentiable simulations emerge and boost as a powerful tool to accommodate various optimization applications crossing graphics and robotics. A typical example is DiffTaichi [42), which created a differentiable programming environment to compute the gradients of physics simulations. A variety of physics simulation algorithms stemming from graphical applications have been adapted to a differentiable framework to facilitate inverse design applications, including fluids [43 44 45 position-basel dynamics 7 cloth [46] 47). deformable objects [48). articulated bodies [49]. object control [50], and solid-fluid coupling systems [39 [1]. While [51]
proposed a method differentiate Lagrangian fluid simulation, optimization of rigid geometry not discussed. Many applications across graphics and robotics have been explored, such as soft-body design and locomotion [52 53 and fluid manipulation [5). However, none of these approaches focused on enabling the inverse design of fluidic device systems in dynamic Navier-Stoke flow.
Computational Design The last decade has witnessed an increasing interest in the design of computational tools and algorithms targeting the digital fabrication of physical systems. A broad
range of applications have been addressed, including the mechanical characters [54 55 56 inflatable
thin shells (57). foldable structures [58] (9). Voronoi structures
and puzzles 661 spinning
objects (62, buoyancy (63] gliders (64) multicopters (65) hydraulic walkers (66. origami robots
1671. articulated robots (68] and multi-materii jumpers 169]. to name just few. Among these applications, the problem of optimizing the shape and control of a 3D printable object to manifest specific mechanical properties and functionalities has drawn particular attention. Examples of designing mechanical properties by optimizing materials include optics 770 (1], mechanical stability
[72,, strength 773. rest shape (74), and desired deformation [75] 76).
6 Conclusions, Limitation and Future Work
In this paper, we proposed fully differentiable pipeline for neural fluidic system control and design, addressing the challenges of complex geometry representation, differentiable fluid simulation, and co-design optimization processes. Our pipeline features liw-dimessional parametric geometry representation and differentiable Navier-Stokes simulator for predicting fluid behavior. We demon- strate the effectiveness of our pipeline a number of complex control design tasks, ranging from different fluidic functional controls to complex neural heart control.
There are certain limitations and avenues for future work. First, the current pipeline assumes the standard Navier-Stokes model, which limits its applicability to Newtonian flow. Extending the framework to handle non-Newtonian flows or multi-physics interactions would be an interesting direction for future research. Additionally, the pipeline relies on parametric representation, which may encounter challenges in navigating complex design spaces with nigh-dimensional or discontinuous parameterizations such as coupling control design with topology optimization. Exploring alternative optimization algorithms or incorporating surrogate models could enhance the efficiency and robust- ness of the optimization process. Furthermore, while we demonstrate the effectiveness of our pipeline in several control and design tasks, additional validation and bench-marking against real-world physical experiments would be valuable to establish the pipeline's 's reliability and generalizability.
10

---

References
[1] G. D. Ntouni, A. Lioumpas, and K. S. Nikita, "Reliable and energy- efficient communications for wireless biomedical implant systems," IEEE Journal of Biomedical and Health Informatics, vol. 18, no. 6, pp. 1848-1856, 2014.
[2] D. Erickson and D. Li, "Integrated microfluidic devices," Analytica Chimica Acta, vol. 507,
pp. 11-26, 2004.
[3] K. Nonami, F. Kendoul, S. Suzuki, W. Wang, and D. Nakazawa, Autonomous flying robols: unmanned aerial vehicles and micro aerial vehicles. Springer Science & Business Media, 2010.
[4] G. Brockman, V. Cheung, L. Pettersson, Schneider, J. Schulman, Tang, and W. Zaremba, "Openai gym," arXiv preprint arXiv: 2016.
[5] Z. Xian, B. Zhu, Z. Xu, H.-Y. Tung, A. Torralba, K. Fragkiadaki, and C. Gan, "Fluidlab: A differentiable environment for benchmarking complex fluid manipulation." ICLR, 2023.
[6] Ma, Du, 7. Zhang, K. Wu, A. Spielberg, R. K. Katzschmann, and W. Matusik, "Diffaqua: A differentiable computational design pipeline for soft underwater swimmers with shape interpolation," " ACM Transactions on Graphics (TOG), vol. 40, no. 4, p. 132, 2021.
[7] T. Du, K. Wu, P. Ma, S. Wah, A. Spielberg, D. Rus, and W. Matusik, "Diffpd: Differentiable projective dynamics," ACM Trans. Graph., vol. 41, no. 2, nov 2021. [Online]. Available: https://dii.rg// 0rt.o//o.orrg0001113344000000000
[8] Y. Li, T. Du, S. Grama Srinivasan, K. Wu, B. Zhu, E. Sifakis, and W. Matusik, "Fluidic topology optimization with anisotropic mixture model." ACM Trans. Graph., nov 2022. [Online]. Available: https://doi.org// http::/000oong111114444444444445 145/3550454.3555522
[9] J. Stam, "Stable fluids," Proceedings of the 26th Annual Conference on Computer Graphics and Interactive Techniques, p. 121-128, 1999.
[10] R. Bridson, Fluid simulation for computer graphics, 2nd ed. Boca Raton, FL, USA: AK
Peters/CRC Press, 2015.
[11] F. H. Harlow and J. E. Welch, "Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface," Physics of Fluids, vol. 8, pp. 2182-2189, 1965.
[12] Y. T. Ng, C. Min, and Gibou, "An efficient fluid-solid coupling algorithm for single-phase flows," Journal of Computational Physics, vol. 228, no. 23, pp. 8807-8829, 2009.
[13] W. Jakob, J. Rhinelander, and D. Moldovan, "pybindll - seamless operability between c++11 and python," " 2016, https:/jithbb.oompyiiid/pipiinnt
[14] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga al., "Py Torch: An imperative style, regh-performance deep learning library," in Neural Information Processing Systems, 2019.
[15] J. Schulman, Wolski, P. Dhariwal, A. Radford, and O. Klimov, "Proximal policy optimization algorithms," 2017.
[16] N. Hansen and S. Kern, "Evaluating the ema evolution strategy multimodal test functions," vol. 3242, 2004, pp. 282-291.
[17] A. McAdams, E. Sifakis, and J. Teran, "A parallel multigrid poisson solver for fluids simulation on large grids," ser. SCA '10, 2010, p. 65-74.
[18] Borrvall and Petersson, "Topology optimization of fluids in stokes flow," International journal for numerical methods in fluids, vol. 41, no. 1, pp. 77-107, 2003.
[19] J. Alexandersen and C. S. Andreasen, "A review of topology optimisation for fluid-based problems." Fluids, vol. 5, no. 1, 2020. [Online]. Available: /www.mapp..cm/23111555511111
https

---

[20] J. K. Guest and J. H. Prévost, "Topology optimization of creeping fluid flows using darcy- stokes finite element," International Journal for Numerical Methods in Engineering, vol. 66, no. 3, pp. 461 1844, 2006.
[21] N. Aage, T. H. Poulsen, A. Sersborg-Hansere and O. Sigmund, "Topology optimization of large scale stokes flow problems," Structural and Multidisciplinary Optimization, vol. 35, no. 2, pp. 175-180, 2008.
[22] V. J. Challis and J. K. Guest, "Level set topology optimization of fluids in stokes flow." In- ternational Journal for Numerical Methods in Engineering, vol. 79, no. 10, pp. 1284-1308, 2009.
[23] Zhou and Q. Li, "A variational level set method for the topology optimization of steady-state navier-stokes flow," Journal of Computational Physics, vol. 227, no. 24, pp. 10 178-10 10 195, 2008.
[24] A. Evgrafov, "Topology optimization of slightly compressible fluids," 7AMM-Journal of Applied Mathematics and Mechanics//eetschriff für Angewandte Mathematik und Mechanik: Applied Mathematics and Mechanics, vol. 86, no. 1, pp. 46 62, 2006.
[25] Y. Deng. Z. Liu, and Y. wu, "Topology optimization of steady and unsteady incompressible navier-stokes flows driven by body forces," Structural and Multidisciplinary Optimization, vol. 47, 2012.
[26] A. Beessoog-HHansen O. Sigmund, and R. B. Haber, "Topology optimization of channel flow problems," Structural and Multidisciplinary Optimization, vol. 30, no. 3, pp. 181-192, 2005.
[27] C. Othmer, E. de Villiers, and H. Weller, "Implementation of continuous adjoint for topol- ogy optimization of ducted flows," in 18th AIAA Computational Fluid Dynamics Conference. Reston, VA, USA: the American Institute of Aeronautics and Astronautics, 2007, p. 3947.
[28] E. Kontoleontos, F. Papoutsis-Kiachagaas A. Zymaris, D. Papadimitriou, and K. Giannakoglou, "Adjoint-based constrained topology optimization for viscous flows, including heat transfer," Engineering Optimization, vol. 45, no. 8, pp. 941-961, 2013.
[29] G. H. Yoon, "Topology optimization for stationary fluid-structure interaction problems using a new monolithic formulation." " International Journal for Numerical Methods in Engineering. vol. 82, no. 5, pp. 591-616, 2010.
[30] W. J. P. Casas and R. Pavanello, "Optimization of Aluid-structure systems by eigenvalues gap separation with sensitivity analysis," Applied Mathematical Modelling, vol. 42, pp. 269-289, 2017.
[31] C. S. Andreasen and O. Sigmund, "Topology optimization of sttructure-itteraeeee prob- lems in poroelasticity," Computer Methods Applied Mechanics and Engineering, vol. 258, pp. 55-62, 2013.
[32] T. Matsumori, T. Kondoh, A. Kawamoto, and Nomura, "Topology optimization for fluid- thermal interaction problems under constant input power," Structural and Multidisciplinary Optimization, vol. 47, no. 4, pp. 571-581, 2013.
[33] K. Yaji. T. Yamada, S. Kubo, K. Izui, and S. Nishiwaki, "A topology optimization method for coupled thermal-fluid problem using level set boundary expressions," International Journal of Heat and Mass Transfer, vol. 81, pp. 878-888, 2015.
[34] C. S. Andreasen, A. R. Gersborg, and O. Sigmund, "Topology optimization of microfluidic mixers," International Journal for Numerical Methods in Fluids, vol. 61, no. 5, pp. 498-513 2009.
[35] M. Mangano, S. He, Y. Liao, D.-G. Caprace, A. Ning, and J. R. R. A. Martins, "Aeroelastic tailoring of wind turbine rotors using high-fidelity multidisciplinary design optimization," Wind Energy Science, Jan. 2023, (in review).
12

---

htt.://o...rr1111114111111706660 [36] Y. Yu, Z. Lyu, Z. Xu, and J. R. R. A. Martins, "On the influence of optimization algorithm and starting design on wing aerodynamic shape optimization." Aerospace Science and Technology. vol. 75, pp. 183-199, Apr. 2018.
[37] A. Jameson, "Optimum aerodynamic design using cfd and control theory," CFD Review, vol. 3, 06 1995.
[38] K. Maute and M. Allen, "Conceptual design of aeroelastic structures by topology optimization." Structural and Multidisciplinary Optimization, vol. 27, no. 1-2, pp. 27-42, 2004.
[39] T. Takahashi, J. Liang, Y.-L. Qiao, and M. C. Lin, "Differentiable fluids with solid coupling for learning and control," AAAI, 2021.
[40] T. Du, K. Wu, A. Spielberg, W. Matusik, B. Zhu, and E. Sifakis, "Functional optimization of fluidic devices with differentiable stokes flow," ACM Trans. Graph., vol. 39, no. 6, Dec. 2020. [Online]. Available: httis:///oii.or org/ I0. 145/3/4444444334177999
[41] A. McNamara, A. Treuille, Z. Popović, and J. Stam, "Fluid control using the adjoint method," ACM Trans. Graph., vol. 23, no. 3, 449-456, aug 2004. [Online]. Available:
[42] Y. Hu, L. Anderson, T.-M. Li, Q Sun, N. Carr, J. Kegan-Kelley, and F, Durand, Taihiiaichi: Differentiable programming for physical simulation," in ICLR, 2020.
[43] Holl, V. Koltun, K. Um, and N. Thuerey, "phiflow: differentiable pde solving framework for deep learning via physical simulations," in Advances in Neural Information Processing Systems (NeurIPS) Workshop, 2022.
[44] P. Holl, V. Koltun, and N. Thuerey, "Learning to control pdes with differentiable physics," 2020.
[45] B. List, L.-W. Chen, and N. Thuerey, "Learned turbulence modelling with differentiable fluid solvers: physics-based loss functions and optimisation horizons," Journal of Fluid Mechanics, vol. 949, Sep. 2022. [Online]. Available: http://dxx.ooi.org// http:/dd.ddo..rr///001////m2.2..
[46] Y. Li, T. Du, K. Wu, Xu, and W. Matusik, "Diffcloth: Differentiable cloth simulation with dry frictional contact," ACM Trans. Graph., mar 2022. [Online). Available: https://doi.or org//0.. 10..1115500 145/3527660)
[47] Y. Li. H.-y Chen, E. Larionov, N. Sarafianos, W. Matusik, and T. Stuyck, "DiffAvatar: Simulation-readd garment optimization with differentiable simulation." in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR). Los Alamitos, CA, USA: IEEE Computer Society, June 2024. [Online]. Available: 1000//..0.00000033550302...000..
[48] Y. Qiao, J. Liang, V. Koltun, and M. Lin, "Differentiable simulation of soft multi-body systems," Advances in Neural Information Processing Systems, vol. 34, 2021.
[49] Y.-L. Qiao, J. Liang. V. Koltun, and M. C. Lin, "Efficient differentiable simulation of articulated bodies," " in International Conference on Machine Learning. PMLR, 2021, pp. 8661-8671.
[50]
"Scalable differentiable physics for learning and control," in ICMI., 2020.
[51] Z. Li, Q. Xu, X. Ye, B. Ren, and L. Liu, "Difffr: Differentiable sph-based fluid-rigid coupling for rigid body control," ACM Trans. Graph., vol. 42, no. 6, Dec 2023. [Online]. Available: https://dbi.orr/10111555//3
[52] Y. Hu, J. Liu, A. Spielberg. J. B. Tenenbaum, W. T. Freeman, J. Wu, D. Rus, and W. Matusik, "Chainqueen: A real-time differentiable physical simulator for soft robotics," 2018. [Online]. Available: https://arxxv..or//a//////1 1100001054
[53] A. Spielberg, A. Zhao, Y. Hu, T. Du, W. Matusik, and D. Rus, TLeearnngginn-th--oon optimiza- tion: End-to-end control and co-design of soft robots through learned deep latent representations," in Neural Information Processing Systems, 2019.
13

---

[54] S. Coros, B. Thomaszewski. G. Noris, S. Sueda, M. Forberg, R. W. Sumner, W. Matusik, and B. Bickel, "Computational design of mechanical characters, ACM Transactions on Graphics (TOG), vol. 32, no. 4, p. 83, 2013.
[55] B. Thomaszewski, S. Coros, D. Gauge, V. Megaro, E. Grinspun, and M. Gross, "Computational design of linkage-based characters," ACM Transactions on Graphics (TOG), vol. 33, no. 4, p. 64, 2014.
[56] D. Ceylan, W. Li, N. J. Mitra, M. Agrawala, and M. Pauly, "Designing and fabricating mechani- cal automata from mocap sequences," ACM Transactions on Graphics (TOG), vol. 32, no. p. 186, 2013.
[57] M. Skouras, B. Thomaszewski. P. Kaufmann, A. Garg, B. Bickel, E. Grinspun, and M. Gross, "Designing inflatable structures," ACM Transactions on Graphics (TOG), vol. 33, no. 4, p. 63, 2014.
[58] Felton, M. Tolley, E. Demaine, D. Rus, and R. Wood, "A method for building self-folding machines," Science, vol. 345, no. 6197, pp. 644-646, 2014.
[59] C. Sung and D. Rus, "Foldable joints for foldable robots," Journal of Mechanisms and Robotics, vol. 7, no. 2, 021012, 2015.
[60] L. Lu, A. Sharf, H. Zhao, Y. Wei, Q. Fan, X. Chen, Y. Savoye, C. Tu, D. Cohen-Or, and B. Chen, "Build-to-last" strength to weight 3d printed objects." ACM Transactions on Graphics (TOG), vol. 33, no. 4, P. 97, 2014.
[61] T. Sun, C. Zheng, Y. Zhang, C. Yin, C. Zheng, K. Zhou, Y. Yue, B. Smith, C. Baatt, Z. Xu el al., "Computational design of twisty joints and puzzles." ACM Trans. Gruph., vol. 34, no. pp. 101-1, 2015.
[62] M. Bächer, E. Whiting, B. Bickel, and O. Sorkine- Hornung, "Spin-it: optimizing moment of inertia for spinnable objects," ACM Transactions on Graphics (TOG), vol. 33, no. 4, p. 96, 2014.
[63] L. Wang and E. Whiting, "Buoyancy optimization for computational fabrication," Computer Graphics Forum (Proceedings of Eurographics), vol. 35, no. 2, 2016.
[64] T. Martin, N. Umetani, and B. Bickel, "Omniad: data-driver omni-directional aerodynamics," ACM Transactions on Graphics (TOG), vol. 34, no. 4. p. 113, 2015.
[65] T. Du, A. Schulz, B. Zhu, B. Bickel, and W. Matusik, "Computational multicopter design," 2016.
[66] R. MacCurdy, R. Katzschmann, Y. Kim, and D. Rus, "Printable hydraulics: A method for fabricating robots by 3d co-printing solids and liquids." 2016.
[67] A. Schulz, C. Sung, A. Spielberg, W. Zhao, R. Cheng, E. Grinspun, D. Rus, and W. Matusik, "Interactive robogami: An end-to-end system for design of robots with ground locomotion," The International Journal of Robotics Research, vol. 36, no. 10, pp. 1131-1147, 2017.
[68] A. Spielberg, B. Araki, C. R. Sung, R. Tedrake, and D. Rus, "Functional co-optimization of articulated robots," in ICRA. IEEE, 2017, pp. 5035-5042.
[69] D. Chen, D. I. Levin, W. Matusik, and D. M. Kaufman, "Dynamics-aware numerical coarsening for fabrication design," ACM Trans. Graph., vol. 34, no. 2017.
[70] M. Hašan, M. Fuchs, W. Matusik, H. Pfister, and S. Rusinkiewicz, "Physical reproduction of materials with specified subsurface scattering," ACM Trans. Graph., vol. 29, no. 4, 2010.
[71] Y. Dong, J. Wang, Pellacini, X. Tong, and B. Guo, "Fabricating spatially-varying subsurface scattering." ACM Trans. Graph., vol. 29, 4, 2010.
[72] O. Stava, Vanek, B. Benes, N. Carr, and R. Mēch, "Stress relief: improving structural strength of 3d printable objects," ACM Trans. Graph., vol. 31, no. 4, 2012.
14

---

[73] Q. Zhou, J. Panetta, and D. Zorin, "Worst-case structural analysis," ACM Trans. Graph., vol. 32, no. 4, 2013.
[74] X. Chen, C. Zheng, W. Xu, and K. Zhou, "An asymptotic numerical method for inverse elastic shape design," ACM Transactions Graphics (Proceedings of SIGGRAPH 2014), vol. 33, no. 4, Aug. 2014.
[75] B. Bickel, M. Bächer, M. A. Otaduy, H. R. Lee, H. Pfister, M. Gross, and W. Matusik, "Design and fabrication of materials with desired deformation behavior," ACM Trans. Graph. vol. 29, no. 4, 2010.
[76] B. Zhu, M. Skouras, D. Chen, and W. Matusik, "Two-scale topology optimization with mi- crostructures," ACM Trans. Graph., vol. 36, no. 4, jul 2017.
A Geometry Representation Implementation
2D Geometry We define a closed 2D surface using connected cubic Bezier curves parameterized in the polar coordinate frame. A point e is defined on the surface to establish the center of a polar coordinate frame. Each Bezier curve spans an are of 2n radians on the polar coordinate plane, N
and its control points are symmetrically placed, each at an angular displacement of 2x. 3N radians from their corresponding curve endpoint. The shape of each Bezier curve, and consequently the overall surface, is manipulated via two scalar parameters, P1 and P2. dictating the polar coordinate distance of the
control points. The i-th cubic Bezier curve is defined by control points Po = (p) cos 0, p1 sinf) + e and = (p} cos 0, phsin0) + c, where e is the reference center point. The endpoints Co. ej computed by ensuring ef e
(1%1 and colinearity of each pair of
(6+1) This representation offers a compact way of defining diverse geometries. pÍ ef. Po
3D Geometry We parameterize a closed 3D surface using 2D surfaces defining the key cross- sections of the geometry along the z-axis of the local object frame, where each 2D surface is parameterized by N cubic Bezier curves. The parametrization includes 2 and ø1, which determines the plane of the first and last cross-section and the parameters for each key cross-section. The i th key cross-section is defined by the center c' and control point parameters p p for each of E [1, 2,....... N). The key cross- sections are assumed to evenly spaced along the z-axis. Then, given 2g I the cross-section of the closed surface at Z is defined by interpolating the centers and control points of all key cross-sections using the interpolation scheme
20
C₂
(10a)
21 20
n
2 20 y
Pj
p
(10b)
21 =0
B Temporal Discretization of the Governing Equation
We build the fluid simulator by leveraging the pperator-splitting method 9 10] Each simulation step comprises of advection, viscosity, and projection.
Advection We employ the semi-Lagrangian advection scheme specified in (Eqs. 11 and 12 to propagate velocity through the fluid domain:
1+12
Vu"
(11)
At/2
x+
tt -1/11/2 V
(12)
At
15

---

Viscosity For incompressible fluid with constant viscosity coefficient, the viscous force density equivalent the product of the Laplacian of velocity and the viscosity coefficient. We employ explicit time integration to update the fluid velocity in response to the viscous force.
m++ n++
v u u+
(13)
A
Projection The projection step involves an update of the pressure and the velocity field (Eq. 14a to ensure the satisfaction of the incompressibility condition (Eq. 14b).
u a=+1__
an+1
-y n+1
(14a)
Al
1 O.
(14b)
On boundaries, the pressure is regulated by two conditions: the Dirichlet boundary condition (Eq. 15a and the non-penetrating Neumann boundary condition (Eq. 15b given computed geometry
velocity u
+1 putl =0, a E 022+ u= n u n. x: E 052
(15a)
(15b)
The pressure field for the subsequent time step p determined by solving the resultant Poisson equation (Eq. I6P.
Al v2
n+l
á m+1
(16)
C Additional Optimization Task Details and Visualization
C.1
Switch
We visualize the initial and optimized design for the amplifier task in Fig. 6
(8....8 Initial Design
Segel
Optimized Design
Figure 6: Fluidic Switch. The switch rotates dynamically across time (dotted lines). The shape of the switch is parameterized as a 2D Polar Bezier, whose parameters, along with the parameters of the rotation signal are subject to optimization. The top and bottom the illustration visualize information from the initial and optimized iteration respectively. For each iteration, we visualize the design geometry (left) and corresponding streamlines of the flow field (right) key-frames evenly sampled across time. We additionally plot the target (green) and outlet velocity norm profile (orange) across time and visualize their difference in grey shaded area.
D Experiment on Fluid Solver Validation Kármán Vortex Street
To further validate the performance of our solver, we conducted an additional experiment simulating the formation of classic Kármán Vortex Street. This experiment was executed at resolution of 512 1024 and illustrates the capability of our solver in capturing complex fluid dynamics phenomena. As shown in Fig. 7 simulate horizontal flow passing around a cylindrical obstacle
16

---

at three distinct kinematic viscosity values: inviscid (v 0), moderate viscosity (2 0.01), and high viscosity (v - 0.1).
Each viscosity setting demonstrates the characteristic vortex shedding pattern associated with Kármán Vortex Street formation. These results confirm solver's ability to replicate this well-known phenomenon and offer insights into the effect of varying viscosity on vortex behavior. This validation experiment supports the accuracy and versatility of the solver across different fluid conditions.
verticit
No Viscosity
Low Viscosity
== 0.002
High Viscosity
0.02
Figure 7: Solver Validation. Visualization of Karman Vortex Street under different viscosity con- ditions. Here, we illustrate the results of the classic Karman Vortex Street test for three different kinematic viscosity values (From top to down our differentiable simulator within a domain size 512 x 1024. Each figure visualizes the vortex patterns, and the results demonstrate how increased viscosity leads to notable change in vortex formation and dissipation.
0.0, 0.002 and 0.02), simulated using
E
Gradient Stability and Solver Steps Statisties
Gradient stability in differentiable physics is well-known challenge, particularly given the potential for gradient explosion vanishing when gradients are accumulated across numerous solver steps. In our approach, however, we have not encountered significant issues with gradient stability. This stability likely due to the accuracy of gradients produced by our framework and the robustness of our numerical solver. To illustrate this, we provide gradient norm statistics over the full course of optimization for three tasks of varying complexity in Table 4 demonstrating consistent gradient magnitudes without evidence of explosion or vanishing. For additional robustness, our implementation includes gradient clipping with threshold 1.0, which can mitigate gradient explosion in particularly challenging scenarios. This technique ensures gradients remain within manageable limits and contributes to the overall stability of our optimization pipeline.
Furthermore, the actual number of solver steps required to advance between frames in our solver depends on the Frirant-Friidricss--Lwwwy (CFL) condition, which is maintained to ensure numerical stability. Additional statistics solver steps over one optimization cycle across various tasks are provided in the upper portion of Table offering further insight into the computational demands and stability characteristies of our approach.
17

---

Metric
Statistic
Shape Identifier
Heart 3D
Gate 3D
Steps
Mean, Std
(11, 0)
(54, 0)
(58, 0)
Min, Max
(0,46, 252.37) (28.04, 38.00)
(1.61, 820.82)
(5.13, (69.21)
Gradient Norm
Mean, Std
(65.96 109.40)
(40.71, 34.36)
Table 4: Statistics the gradient norm and step count over the full course of optimization for three tasks of varying complexity.
Table 5: Gradient Validation for Shape Identifier Task
Gradient Values by Parameter (P1 to P6)
P1
P2
P3
P4
P5
P6
Analytic
0.048
0.394
0.314
0,046
0.133
0.067
Finite Diff
0.048
0.396
0.314
0.046
0.134
0.066
Abs Diff
-1.7e-4
-2.1e-3
5.2e-4
-3.4e-4
-1.4e-3
7.2e-4
Elem Err
0.004
0.005
0.002
0.007
0.011
0.011
Gradient Values by Parameter (P7 to P11)
P7
P8
P9
P10
P11
Analytic
0.087
0.008
-0.022
0.540
-0.058
Finite Diff
0.086
0.008
-0.025
0.540
-0.059
Abs Diff'
2.1e-4
-9.6e-5
2.2e-3
2.4e-4
9.5e-4
Elem Err
0.003
0.012
0.098
0.000
0.016
F Experiment on Gradient Validation
To ensure the correctness of the gradients in our differentiable simulation framework, we validated the analytical gradients all kernels, functions, and the entire simulation and optimization pipeline using finite difference approximations, Specifically, we employed the central difference method with step size of 1.2 10 to approximate the gradients numerically and compared them with the analytical gradients calculated by our solver.
In this validation experiment, consider the end-to-end gradients for the Shape Identifier Task, which involves optimizing over 1l parameters. The analytical gradients, finite difference gradients, their absolute differences, and element-wise errors are reported in Table For this task, observed a relative error of 0.0047 for the gradient vector, confirming the high accuracy of our analytical gradients.
18

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
The claim reflects paper's contribution and scope accurately.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results be expected to generalize other settings.
is fine include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
We discussed the limitations of the work in Sec. 6
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification. asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low images are taken in low lighting. Or speech-to-text system might be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question:
For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [NA]
19

---

Justification:
We did not include theoretical results.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer:[Yes
Justification:
We discussed the experimental setup in details.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make clear how to reproduce that algorithm. (b) the contribution is primarily a new model architecture, the paper should describe the architccture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. the case of closed-source models, may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
20

---

Answer:[Yes
Justification:
We are releasing the experiment data and code on github.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.cm
ooiiioooooyyylloool for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
/ /pipps.cc/pulli///uuii wwwiii/gguie///eeeeeeoeic CoodeSubeiiieeee ssi wswonnolli onnolicy)))) for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification:
We provide experimental setup in details.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: We provide loss-iteration curves to report experiment details.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, at least for the experiments that support main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the error bar is the standard deviation the standard error of the mean.
21

---

It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of enrors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: Yes]
Justification: We discussed the experimental setup in the paper.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g.. preliminary or failed experiments that didn't make into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https:: https//cccccci.i/tticcccccceeeee
Answer: [Yes]
Justification:
The research conform NeurlPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is a special consid- eration due to laws or regulations their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: There is societal impact of the work performed.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to
22

---

generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
there are negative socictal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: The paper poses such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [NA]
Justification:
The paper does not use existing assets.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, wwppppswipteeeee whw.ree.wthhhe....c.. com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
this information is not available online, the authors are encouraged to reach out to the asset's is creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
23

---

Answer: [NA]
Justification: Our paper does not release new assets.
Guidelines:
The answer NA means that the paper does release new assets.
Researchers should communicate the details the dataset/coddemoodll as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
24


Reviews:
[{"review_id": "3vIP1iErSp", "reviewer": "Reviewer_grLC", "summary": "The paper aims at a fully automated pipeline for devising neural controls for complex fluidic systems with dynamic boundaries. The system consist of externally driven soft boundaries and internal complex flow behaviors.\n\nThe proposed framework contains a differentiable geometry representation, a differentiable fluid simulator with solid-fluid interface handling, and a control-shape co-design algorithm using gradient-based optimization.\n-  The 3D surface is represented by a list of 2D surfaces (represented by a set of connected Bezier curves of their control points) defining the first and last cross-section, and each key cross-sections of the geometry. More complex 3D geometries are represented by union and intersection of sub-geometries.\n- The fluid simulator leverages the operator-splitting method. A single simulation step comprises three sub-steps: advection, viscosity, and projection. It uses MAC grid. \n- The back-propagation process is used to construct the gradients of the geometry surface w.r.t. the loss function.\n\nThe method is evaluated on multiple tasks, from simpler 2D tasks to complex 3D tasks (such as, \nneural gate controller and artificial heart design).", "questions": "- It might be useful to demonstrate whether gradients explode or vanish during the iterative back-propagation process.\n- In the experiment studies the effect of random initialization in our fluid optimization tasks, the Shape Identifier task (with #Parameter = 10) is used. Does the same conclusion hold for more complex tasks with a greater number of parameters?\n- How sensitive are the design results to the grid resolution? Given that the grid used for generating the design is relatively small for fluid simulations, can the designed system achieve the same level of accuracy when evaluated on a higher resolution grid?\n- How does the method perform when using a higher resolution grid (but not so high that the simulator cannot finish within an acceptable time)?\n- In Table 2: Time Performance, what's the time unit?", "limitations": "The author have discussed the limitations in the paper.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The paper contains a fair amount of work to develop a differentiable fluid simulation and optimization framework, which includes a differentiable geometry representation, a differentiable fluid simulator with solid-fluid interface handling, and a gradient-based optimization framework.\n- The proposed methods are evaluated on multiple cases from simple 2D examples to complex 3D tasks. Besides, it's nice to see the performance profiling and comparison results in the paper.", "weaknesses": "It would be nice to see the framework open-sourced.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the constructive suggestions. Below we provide responses to individual questions.\n\n**1. It would be nice to see the framework open-sourced.**\n\nWe will open-source our code upon acceptance. We have provided the anonymized version of our code at https://anonymous.4open.science/r/DiffFluidRebuttal-1D24/.\n\n**2. It might be useful to demonstrate whether gradients explode or vanish during the iterative back-propagation process.**\n\n * Thank you for your insightful point. The issue of gradient explosion or vanishing is indeed a common problem in differentiable physics. However, we have not observed significant problems in our case, likely due to the clean gradients we provide and the stability of our numerical solver. \n * We have included additional statistics (min,max,mean, standard deviation) on the gradient norm over the full course of optimization for three tasks with varying complexities in Table 1. Our results do not show evidence of gradient explosion or vanishing.\n* In practice, in case a gradient explosion is encountered, we implement gradient clipping of the gradient of the network to 1.0. We will include this detail in the implementation section of our paper to enhance the clarity and robustness of our approach.\n\n\n**3. In the experiment studies the effect of random initialization in your fluid optimization tasks, the Shape Identifier task (with #Parameter = 10) is used. Does the same conclusion hold for more complex tasks with a greater number of parameters?**\n\nYes, the same conclusion holds for all of our tasks. In Fig 2 we present an additional verification with the 3D heart controller task which has 7.1k neural network parameters. We initialize the network parameters with 5 different random seeds, and observe the same convergence across seeds even with different initial behaviors, showing the robustness of our method.\n\n**4. How sensitive are the design results to the grid resolution? Given that the grid used for generating the design is relatively small for fluid simulations, can the designed system achieve the same level of accuracy when evaluated on a higher resolution grid? How does the method perform when using a higher resolution grid (but not so high that the simulator cannot finish within an acceptable time)?**\n\n* Thank you for raising this important question. In pure forward simulation, higher grid resolutions are necessary to achieve greater accuracy, particularly at the solid-fluid interface, which directly impacts optimization results. This is because the same continuous geometry boundary can yield different simulation results depending on the resolution, the objective values—and consequently the optimized designs—will vary across different resolutions. \n* To illustrate this, we have conducted an additional experiment to analyze the effect of grid resolution on design outcomes, where we first optimize for the Shape Identifier task under 128x128, then take the optimized values to 512x512, which evaluates to a loss value of 3.11. In comparison, a randomly initialized design under 512x512 evaluates to a loss value of 245.1 and optimizes to a design with loss 0.968.  The results demonstrate that while higher resolutions generally improve accuracy, they also introduce variations in the optimal design due to the increased detail captured at the boundary. \n\n\n**5. In Table 2: Time Performance, what's the time unit?**\n\nThe time reported is in Seconds (s)", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "WIidwvNFvr", "reviewer": "Reviewer_ZAVx", "summary": "The paper proposes a new set of utilities for experimenting with system design and control of viscous fluid flows on deformable domains. The contributions include (a) a differentiable NSE solver, (b) Bezier curve-based geometry parametrization, (c) an algorithm to jointly optimize a control and design objective, (d) a reinforcement learning environment interface, and (e) benchmark cases with baseline results. The results demonstrate the superiority of the differentiable solver framework over existing gradient-free approaches and pave the road for exciting research, e.g., in medical applications.", "questions": "- Did you consider submitting to the Datasets & Benchmarks track instead of the general conference? Your contribution feels more like a benchmarking suite than machine learning research.\n- In the caption of Table 1, you mention that \"# Frames\" is not the same as the number of solver steps. What is the actual number of solver steps between two frames? Accumulating gradients over hundreds of steps often leads to exploding gradients, so why don't you have these issues if you stick to the CFL number? \n- Related to the previous question, using learned surrogates as PDE solvers has demonstrated great potential in overcoming the exploding gradients issue by simply learning to do as much as 200x larger time integration steps [1]. Does your framework allow for replacing the NSE solver with a learned surrogate, e.g. a U-Net or FNO?\n\n---\n[1] Allen et al., \"Inverse Design for Fluid-Structure Interactions using Graph Network Simulators\", NeurIPS 2022", "limitations": "The authors have fairly assessed the limitations of the proposed approach.", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "- Implementing an efficient PDE solver from scratch in C++ and CUDA and then analytically deriving gradients is a great service to the community.\n- The manuscript is well written and provides a good balance between theory and experimental evidence.\n- Efficiently parametrizing complex geometries has a great potential for bridging the gap between toy experiments and real-world problems.", "weaknesses": "- **A. Proper figures**: make the figures on pages 3 and 4 proper figures with captions. Currently, the only way to refer to them is as \"the figure on paper 3 ...\"\n- **B. Solver validation**: having some experience implementing numerical solvers, I know that validating the implementation is crucial. Could you add an appendix on validating the solver on an established fluid mechanics benchmark, like channel flow? The same applies to the gradients, which are typically validated by comparing them against a finite difference numerical approximation.\n- **C. Minor issues**: Line 131 \"can derived\"?; Line 155 \"of from\"?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for appreciating our work and raising constructive suggestions.\n\n**1. Proper figures: make the figures on pages 3 and 4 proper figures with captions\"**\n We will make both inset figures proper figures in the revision.\n\n**2. Solver validation: having some experience implementing numerical solvers, I know that validating the implementation is crucial. Could you add an appendix on validating the solver on an established fluid mechanics benchmark, like channel flow?**\n\nThank you for your suggestion. We have added an experiment validating our solver using the classic Karman Vortex Street under a resolution of 512x1024. In Fig.1, we show the vortex pattern formed by a horizontal flow passing through a cylinder at three different kinematic viscosity values: $\\nu = 0.0$ (inviscid fluid), $\\nu = 0.002$, and $\\nu = 0.02$ to recreate the classic phenomenon. We will add this experiment to the Appendix.\n\n\n**3. The same applies to the gradients, which are typically validated by comparing them against a finite difference numerical approximation.** \n\nIndeed, it is important to validate the correctness of the gradients. During the development of our differentiable simulation, we validated the correctness of our analytical gradients for all kernels, functions, and the simulation and optimization pipeline using finite difference approximations. Below, we include the validation experiment for the end-to-end gradients of the Shape Identifier Task, which has a total of 11 optimization parameters. We compute the finite difference gradient using the central difference with a step size of 1.2e-5. Below we report the two gradients, their element-wise, as well as the vector error and difference.\n\n|||||||||||||\n|---|---|---|---|---|---|---|---|---|---|---|---|\n|Analytical|0.0483|0.3935|0.3142|0.0459|0.1327|0.0670|0.0866|0.0081|-0.0224|0.5404|-0.0579|\n|Finite Difference|0.0484|0.3956|0.3137|0.0462|0.1341|0.0663|0.0864|0.0082|-0.0246|0.5402|-0.0589|\n|Absolute Difference|-1.722e-04|-2.090e-03|5.196e-04|-3.385e-04|-1.415e-03|7.247e-04|2.143e-04|-9.613e-05|2.195e-03|2.352e-04|9.514e-04|\n|Element-wise Error|0.0036|0.0053|0.0017|0.0074|0.0107|0.0108|0.0025|0.0119|0.0981|0.0004|0.0164|\n\n$\\lVert \\text{diff} \\rVert$: 0.0036\nRelative Error: 0.0047\n\n**4. Did you consider submitting to the Datasets & Benchmarks track instead of the general conference? Your contribution feels more like a benchmarking suite than machine learning research.**\n\nWhile our submission does include a suite of tasks for inverse problems in fluid environments, we believe it goes beyond merely being a benchmarking suite. Our paper introduces a novel framework that provides a comprehensive parametric geometry pipeline for specifying complex geometries and offers gradients through the geometry layer. This combination, along with our efficient solver and gradient computation implementation, enables optimization and manipulation tasks on complex shape boundaries that were previously unachievable. These contributions position our work as an advancement in the differentiable simulation community, particularly in its application to robotics, computational fluid dynamics (CFD), and inverse design. Therefore, we believe the general conference track is more appropriate for our submission as it highlights the innovative aspects of our research.\n\n\n**5. In the caption of Table 1, you mention that \"# Frames\" is not the same as the number of solver steps. What is the actual number of solver steps between two frames? Accumulating gradients over hundreds of steps often leads to exploding gradients, so why don't you have these issues if you stick to the CFL number?**\n* Thank you for your insightful question. The issue of gradient explosion or vanishing is indeed a common problem in differentiable physics. However, we have not observed significant problems in our case, likely due to the clean gradients we provide and the stability of our numerical solver. We provide the statistics of the gradient norm over the full course of optimization on three tasks with varying complexities in Table 1. Our results do not show evidence of gradient explosion or vanishing. In practical engineering scenarios, if a gradient explosion is encountered, we implement gradient clipping and clip the gradient of the network to 1.0. We will include this detail in the implementation section of our paper to enhance the clarity and robustness of our approach.\n* To address your question regarding the number of solver steps, we have provided additional statistics on the actual number of solver steps in the upper part of Table 1 over one full course of optimization on three tasks. \n\n**6. ...using learned surrogates as PDE solvers has demonstrated great potential in overcoming the exploding gradients issue by ... 200x larger time integration steps [1]. Does your framework allow for replacing the NSE solver with a learned surrogate, e.g. a U-Net or FNO?**\n\nOur geometry and optimization framework is designed to be orthogonal to our solver module, allowing for integration with learned NSE solvers. One advantage of our PDE solver over learned surrogates is its ability to solve any given combination of geometry boundary and initial conditions without prior training. In contrast, a learned alternative typically requires extensive training on a large dataset with various boundary conditions and geometries to achieve sufficient generalizability for geometry optimization. The design and control tasks that we highlight in this paper aim to explore novel geometry designs and control in unseen simulation conditions, which a learned surrogate might struggle to do and remains a challenging and active area of research.\nWe acknowledge that learned surrogate models have the potential to enhance traditional differentiable physics simulators by enabling larger time steps, and we view this research direction as complementary to our work and a promising area for future development.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "0ys98X4eIz", "reviewer": "Reviewer_QmuX", "summary": "This paper introduces a novel approach to fluidic system design and control using differentiable simulation. The authors propose a method that leverages gradient-based optimization to enhance performance and accuracy in fluid dynamics applications.", "questions": "The majority of the novelty comes from an efficient CUDA kernel implementation but the specifics are omitted. \nWhat specific design considerations are taken that enable this method to outperform existing solvers?", "limitations": "Although the authors outline the limitations of their method for other fluid models, they fail to address significant limitations encountered in many differentiable fluid simulations. Specifically,during backpropagation, such as the need to efficiently checkpoint the solution or the dependence of the amount of memory needed.", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The paper demonstrates the effectiveness of their approach by comparing it to gradient-free methods such as Proximal Policy Optimization (PPO) and Covariance Matrix Adaptation Evolution Strategy (CMA-ES).\nThe authors effectively highlight the advantages of having access to gradients, allowing for better optimization due to more information about the underlying function. In contrast, gradient-free methods approximate these gradients, which can be less efficient.\nThe approach is particularly beneficial in scenarios where the underlying dynamics are smooth, as opposed to scenarios better suited for gradient-free methods.", "weaknesses": "The concept of backpropagation through time in reverse order has been extensively explored in neural ordinary differential equations (neural ODEs) and various adjoint methods literature.\nThe paper appears to be closely related to existing work but utilizes a different discretization or numerical method. The claim of differentiability may be redundant as the equations themselves, or their adjoint operators, are inherently differentiable.\nThe approach resembles the process of writing custom Jacobian-vector products (JVPs), which may not significantly advance the field.Although the authors outline the limitations of their method for other fluid models, they fail to address significant limitations encountered during backpropagation, such as the need to checkpoint the solution or the increasing memory use as a function of simulation length.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We are grateful for the reviewer’s insightful comments and questions.\n\n**1. The majority of the novelty comes from an efficient CUDA kernel implementation but the specifics are omitted. What specific design considerations are taken that enable this method to outperform existing solvers?**\n* Because differential operators and interpolations access neighboring cells of a spatial point in all directions, we divide the whole simulation domain into cubic blocks, and each cubic block corresponds to a CUDA block. When launching a CUDA kernel, we first load the simulation data of a block into shared memory, and perform calculation efficiently on shared memory. To increase the memory throughput of the data transfer between global memory and shared memory, the simulation data of a cubic block is stored consecutively on global memory. \n* While we adopt a matrix-free MGPCG solver for Poisson, Phiflow uses a CG solver. Therefore, our solver has a much better convergence rate than theirs. To implement an efficient MGPCG solver, we also devised a hierarchical grid data structure on GPUs and customized CUDA kernels for prolongation and restriction operations between coarse and fine grid. Our GPU solver is highly optimized and specially tuned for fluid simulation tasks, which will be open-source.\n\n**2. ...such as the need to efficiently checkpoint the solution or the dependence of the amount of memory needed.**\n\n* We are glad that the reviewer raised this limitation. Gradient checkpointing is a common technique in deep learning applications involving large models. However, our method is specifically designed for differentiable fluid simulation, where all the gradients have been manually derived for efficiency instead of relying on automatic differentiation like DiffTaichi. As a result, our method does not require the automatic differentiation workflow to store all the intermediate activations inside each timestep. In particular, our adjoint derivation of the project solve step is independent of the number of iterations in the solver, which makes our method a way better choice for advection-projection fluid simulation than any existing baselines. \n* To support our claim, we have conducted additional experiments to include a comparison with DiffTaichi, which uses automatic differentiation to derive gradients. We implemented a Conjugate Gradient (CG) iterative solver for the projection step, which is the most time-consuming component of each fluid solve, in DiffTaichi and compared its time and memory performance with our projection solver across four different resolutions in a 3D optimization scenario. The results of this comparison are detailed in Table 2.\n* To summarize, our solver uses significantly less memory than DiffTaichi, with the advantage becoming more pronounced at higher grid resolutions. For instance, at a 64x64x64 resolution, DiffTaichi uses 12 times more memory. This is because DiffTaichi relies on automatic differentiation and needs to store the computational graph and values for all CG iterations, whereas our adjoint gradient method does not require storing these values. This makes our solver more suitable for long-term optimization at high resolutions, where efficient memory usage is crucial. In conclusion, we acknowledge that gradient checkpoint is necessary when the computational graph grows, but our method has optimized the memory consumption to the best extent.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "0mbFqlO64x", "reviewer": "Reviewer_ueL6", "summary": "The paper details a framework for design and control of complex fluidic systems with soft and hard boundaries. The paper specifically develops a fully differentiable pipeline that is able to optimize the design of the mesh to achieve the specified fluid control goal. The proficiency of the proposed framework has been demonstrated on multiple design and fluid control tasks like: (i) amplifying flow of a fluid by changing the shape of an initial boundary (ii) optimizing the position of a \"neural gate\" to satisfy the target outflow requirement (iii) designing an artificial heart (starting from an initial parameterized representation) to achieve a low-error match w.r.t a target simulation. The diversity of the tasks as well as task complexity and good results make the contributions in the paper worth considering. However, the paper lacks crucial details without which it is hard to fully appreciate the effectiveness of the proposed methods, hence making it hard to support publication of the paper in the present state.", "questions": "See questions in `Weaknesses` section.", "limitations": "One improvement that can be made would be to highlight the details of the initial design specifications for the various tasks and enumerating some general prescriptions about generating good initial designs to be supplied to the pipeline.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "* The paper tackles an interesting, challenging and important problem of developing a (fully differentiable) framework for optimal design and control of fluid flows in challenging contexts.\n\n\n* The results presented in the paper have been demonstrated on challenging tasks (i.e., not the usual \"toy\" problems tackled by other such papers) and as presented seem to be effective on the design and control tasks investigated.", "weaknesses": "* Although the results seem interesting and the problem the paper tackles is challenging, the paper lacks articulation of crucial details making the extent of contribution hard to appreciate. Some crucial aspects of the paper need further clarification.\n\n> Q1. How is the \"initial parametric geometry\" (required as input by the framework) generated? What are the feasibility constraints as far as this initial geometry is concerned? Where can a new user of this framework obtain such an initial geometry to supply to the system? An analysis of this or at least an articulation for usability by users is necessary.\n\n> Q2. Further details of the experiment setup are necessary. For example: what is the design of the PPO, CMA-ES models employed? What is the reward function for PPO and what were the architectures of the learnable components therein? \n\n> Q3. Why was `PhiFlow` chosen as the framework for comparison? Wouldn't `DiffTaichi`[1] be a more apples-to-apples comparison for measuring speedup, considering that the C++ implementation of underlying libraries similar to the current proposed framework?  \n\n* Section 2.2  needs further detailing for a machine learning audience who may not be familiar with `Constructive Solid Geometry`.  \n\n* Code pipeline has not been made available for review. Hence, the results are currently not reproducible. It is imperative to make code available during the review stage.\n\n### `References`\n1. Hu Y, Anderson L, Li TM, Sun Q, Carr N, Ragan-Kelley J, Durand F. Difftaichi: Differentiable programming for physical simulation. arXiv preprint arXiv:1910.00935. 2019 Oct 1.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the careful review and insightful questions.  \n\n**1. How is the \"initial parametric geometry\" generated? What are the feasibility constraints...?**\n* Thank you for bringing up this important question. Usability is crucial for allowing users to specifying complex shapes enabled by our framework. We therefore designed the geometry specification process similar to existing CAD software. \n* For instance, a 3D component can be created by extruding a 2D shape along an axis, and users can union or intersect 3D models to form complex structures. Each 2D component is defined using a polar Bézier specification. This approach ensures that users with basic CAD software experience can effectively use our model. \n* For example, the heart model presented in our paper was created by specifying each heart component using six cross-sections, each made of circles with various radii and center positions. This method makes the geometry specification process intuitive and accessible. * By leveraging familiar CAD concepts, we aim to make our framework user-friendly and accessible to researchers and practitioners. We will also provide detailed documentation, visualization tools, and examples to further assist new users in generating the initial parametric geometry needed for their specific applications when releasing the codebase.\n\n**2. Code pipeline has not been made available for review.**\n\nWe will open-source our code upon acceptance. We have provided the anonymized version of our code at https://anonymous.4open.science/r/DiffFluidRebuttal-1D24/ for the review stage.\n\n**3. Sec 2.2 needs further detailing for a ML audience who may not be familiar with CSG.**\n\nThank you for the suggestion. We will add a more detailed description of CSG in our revised manuscript. Briefly speaking, CSG defines a hierarchical tree structure to represent volumetric geometry. Each leaf node in the tree denotes a geometric primitive (e.g., the Bezier surfaces in Sec. 2.2), and each intermediate node up to the root forms a new volumetric geometry by applying boolean operators to its child nodes, e.g., taking the union of the volumetric geometries represented by the child nodes. The root node represents the final geometry, which is formed by recursively applying boolean operators (typically unions, intersections, and subtractions) to a set of geometric primitives stored in the leaf nodes.\n\n**4. Further details of the experiment setup are necessary. What is the design of the PPO, CMA-ES models employed? What is the reward function for PPO and what were the architectures of the learnable components therein?**\n\nOur PPO implementation is based on the open-sourced code of pytorch-a2c-ppo-acktr. Specifically, the actor and critic networks are both two-layer MLPs with tanh activations and a hidden layer of 64 neurons. For both controller tasks, we set the observation to be the same as the input to our gradient-based closed-loop controller. The reward function is the negative comparative of the loss function of the gradient-based closed-loop controller, plus a positive component to avoid negative rewards. Our CMA-ES implementation is based on the open-source Nevergrad. We set the population size to 10 and the metric to be the same as the loss function used for our method. We will include the details of the baseline in the Appendix.\n\n**5. Why was PhiFlow chosen as the framework for comparison? Wouldn't DiffTaichi[1] be a more apples-to-apples comparison for measuring speedup, considering that the C++ implementation of underlying libraries similar to the current proposed framework?**\n\n*Thank you for your valuable feedback. We chose PhiFlow for our comparison because it is an open-source, differentiable simulator specifically optimized for fluid simulation, which aligns closely with our objectives. While DiffTaichi is indeed a powerful differentiable programming tool, it serves a broader range of applications similar to Warp or JAX.\n* However, we recognize the importance of providing a comprehensive comparison for practitioners. We have conducted additional experiments to include a comparison with DiffTaichi. In these experiments, we implemented a Conjugate Gradient (CG) iterative solver for the projection step, which is the most time-consuming component of each fluid solve, in DiffTaichi and compared its time and memory performance with our solver in 3D under 4 resolutions. The results of this comparison are detailed in Table 2.\n* To summarize, our implementation is approximately 3-4 times faster than DiffTaichi in both forward and back propagation. This performance gain is likely due to our efficient CUDA kernel tailored for the Laplacian operator and the rapid convergence of our MGPCG solver. Additionally, our solver uses significantly less memory than DiffTaichi, with the advantage becoming more pronounced at higher grid resolutions. At 64x64x64, DiffTaichi uses 12 times more memory. This is because DiffTaichi relies on automatic differentiation and needs to store the computational graph and values for all CG iterations, whereas our adjoint gradient method does not require storing these values. This makes our solver more suitable for long-term optimization at high resolutions, where efficient memory usage is crucial.\n\n**6. One improvement that can be made would be to highlight the details of the initial design specifications for the various tasks and enumerating some general prescriptions about generating good initial designs...**\n\nThank you for your valuable suggestion. For training fluid controllers, such as in the neural gate and neural heart tasks, we have found that initializing the network to output small actions helps achieve stable initial simulations and facilitates faster convergence of the network. In our implementation, this is accomplished by initializing all network layers with orthogonal matrices scaled to a small norm. We will include additional details on initial design in the Appendix.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We thank all reviewers and the AC for their time and effort in reviewing and for insightful comments to strengthen our work. Besides the responses to individual reviewers, here we would like to highlight our contributions and new quantitative/qualitative results added in the rebuttal.

1. **Contributions**
    1. **[Motivation]** Our manuscript tackles an important, challenging, and significant problem in the optimal design and control of fluid flows in complex contexts [Reviewer ueL6] and effectively bridges the gap between toy experiments and real-world problems [Reviewer ZAVx].
    2. **[Method]** Our method is novel, implementing an efficient PDE solver from scratch in C++ and CUDA with analytically derived gradients, providing a great service to the community [Reviewer ZAVx, Reviewer grLC]. The differentiable fluid simulation and optimization framework includes a differentiable geometry representation and fluid simulator with solid-fluid interface handling [Reviewer grLC].  
    3. **[Experiments]** Our extensive experiments cover a range of scenarios from simple 2D examples to complex,challenging 3D tasks [Reviewer grLC], not just usual "toy" problems [Reviewer ueL6]. We also provide performance profiling and comparison results [Reviewer grLC] and highlights the advantages gradient-based optimizations compared to gradient-free methods [Reviewer QmuX].
    4. **[Presentation]** Our manuscript is well-written, and provides a good balance between theory and experimental evidence [Reviewer ZAVx, Reviewer grLC].


2. **New Results**
    1. **[Code: Open-Source Release]**
    https://anonymous.4open.science/r/DiffFluidRebuttal-1D24/ We have anonymously released our code open-source for the review stage and will officially release it upon acceptance, enabling the community to reproduce our results and build upon our work.
    2. **[Experiment: Solver Validation]**. We validated our fluid solver using the classic Karman Vortex Street test.
    3. **[Experiment: Gradient Validation]**. We presented a gradient validation test to compare our analytical gradient to finite-difference gradient.
    4. **[Experiment: Gradient Norm Analysis]**. We provide statistics on gradient norm on multiple tasks during optimization to study whether gradients diminish/explode.
    5. **[Experiment: Step Number Statistics]**. We provide statistics for solver steps on multiple tasks.
    6. **[Experiment: Effect of Initialization on Complex Tasks]**. We provide an additional experiment on the effect of initialization on the 3D neural controller task with 7k parameters.
    7. **[Experiment: Memory & Time Comparison with DiffTaichi]**. We compare the memory usage and computational time of our solver with DiffTaichi, highlighting the advantages of our implementation in terms of efficiency and scalability.
    8. **[Experiment: Effect of Grid Resolution on Optimization]**. We explore how varying grid resolutions impact the optimization outcomes.
    9. **[Implementation Details: GPU Code Framework & Kernel Acceleration]**. We detailed our GPU code framework and kernel acceleration techniques.

Meta-review: The paper describes a neural network-based framework for the design and control of fluidic systems that involve dynamic (deformable) boundaries. The framework includes a differentiable Navier-Stokes solver, a parametric (Bezier) representation of the geometry, and an algorithm for joint optimization of control and design. The optimization and evaluation utilize a proposed gym-like RL environment along with a series benchmark experimental scenarios. Experiments, including those that involve the proposed benchmarks, reveal the advantages of the framework.

The paper was reviewed by four referees whose overall assessments of the paper all lean positive. Several reviewers find that the problem is interesting and that the algorithmic implementation (i.e., the C++/CUDA PDE solver and the Bezier parameterization) provide a valuable contribution to the community. Others emphasize the thoroughness of the experimental evaluation. There were some concerns about the significance of backpropagation in time, with one reviewer commenting that it has been explored extensively for neural ODE solvers. In this way, the two reviewers find that the contributions relative to existing work is not sufficiently clear (i.e., with regard to the CUDA implementation). The same reviewer emphasizes that the paper does not adequately discuss the limitations encountered during backpropagation with regard to simulating differentiable fluids. In their response to the reviewers, the authors provide more detail regarding their CUDA kernel implementation, which helps to clarify the contribution as the implementation of a custom fluidic solver, as well as elaborate on issues related to backpropagation. However, as the reviewer points out, the author response could have better conveyed the paper's contributions and its relevance to the community. Meanwhile, other reviewers comment on the need to provide a validation of the solver, which the authors do in their rebuttal, and to make the code publicly available, which the authors clarify will be done if and when the paper is published.

Overall, the quality of the paper improved as a result of the initial reviews and the author feedback, as indicated by the more favorable reviewers maintaining their high scores and the two less favorable reviewers increasing their score. The AC agrees with Reviewer QmuX that the paper could be improved further by being more explicit about the contributions relative to existing methods and the relevance to the community.
Decision: Accept (poster)