NeuralFluid: Nueral Fluidic System Design and Control with Differentiable Simulation

We thank the reviewer for the careful review and insightful questions.

1. How is the "initial parametric geometry" generated? What are the feasibility constraints...?

• 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.
• 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.
• 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.

2. Code pipeline has not been made available for review.

We 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.

3. Sec 2.2 needs further detailing for a ML audience who may not be familiar with CSG.

Thank 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.

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?

Our 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.

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?

*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.

• 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.
• 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.

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...

Thank 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.

Dear reviewer ueL6, we wanted to check if there are any remaining questions or concerns about our rebuttal. Since the author-reviewer discussion period ends tomorrow, we’d appreciate any final feedback you might have.

Thank you for your time and attention.

Dear Reviewer,

As the discussion deadline approaches, we are eager to address any remaining concerns you may have. We kindly request you to review our responses. If there is anything further we can do to facilitate the evaluation process, please let us know. We also respectfully ask that you reconsider our score based on the rebuttal provided.

Thank you for your attention.

Best regards,

Authors

We are grateful for the reviewer’s insightful comments and questions.

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?

• 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.
• 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.

2. ...such as the need to efficiently checkpoint the solution or the dependence of the amount of memory needed.

• 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.
• 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.
• 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.

Dear reviewer QmuX, we wanted to check if there are any remaining questions or concerns about our rebuttal. Since the author-reviewer discussion period ends tomorrow, we’d appreciate any final feedback you might have.

Thank you for your time and attention.

Dear Reviewer,

As the discussion deadline approaches, we are eager to address any remaining concerns you may have. We kindly request you to review our responses. If there is anything further we can do to facilitate the evaluation process, please let us know. We also respectfully ask that you reconsider our score based on the rebuttal provided.

Thank you for your attention.

Best regards,

Authors

We thank the reviewer for appreciating our work and raising constructive suggestions.

1. Proper figures: make the figures on pages 3 and 4 proper figures with captions" We will make both inset figures proper figures in the revision.

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?

Thank 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.

3. The same applies to the gradients, which are typically validated by comparing them against a finite difference numerical approximation.

Indeed, 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.

$\lVert \text{diff} \rVert$: 0.0036 Relative Error: 0.0047

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.

While 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.

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?

• 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.
• 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.

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?

Our 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. We 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.

Dear reviewer ZAVx, we wanted to check if there are any remaining questions or concerns about our rebuttal. Since the author-reviewer discussion period ends tomorrow, we’d appreciate any final feedback you might have.

Thank you for your time and attention.

Thank you for your detailed feedback and for providing a thorough suggestion of the validation process in fluid mechanics.

Regarding your points:

1. Numerical Viscosity in Semi-Lagrangian Scheme: We acknowledge that our semi-Lagrangian scheme introduces numerical dissipation [1], which explains the results in Figure 1 in the rebuttal PDF for the chaotic turbulence case. We will add this in the discussion section of our method.

2. Validation with Reference Data: Thank you for the additional feedback. To further validate our solver numerically, we focused on validating the mid-viscosity scenario that we provided in the rebuttal PDF, corresponding to a Reynolds number of 209. For this case, we calculated a Strouhal number of 0.16, compared to the reference value of 0.18 from Irvine 1999. The observed discrepancy can be attributed to the numerical dissipation introduced by our scheme.

We appreciate your insights, which have helped us better articulate the limitations and implications of our approach. Please let us know if you have further questions.

[1] Jos Stam. 1999. Stable fluids. In Proceedings of the 26th annual conference on Computer graphics and interactive techniques (SIGGRAPH '99). ACM Press/Addison-Wesley Publishing Co., USA, 121–128. https://doi.org/10.1145/311535.311548

We thank the reviewer for the constructive suggestions. Below we provide responses to individual questions.

1. It would be nice to see the framework open-sourced.

We will open-source our code upon acceptance. We have provided the anonymized version of our code at https://anonymous.4open.science/r/DiffFluidRebuttal-1D24/.

2. It might be useful to demonstrate whether gradients explode or vanish during the iterative back-propagation process.

• 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.
• 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.
• 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.

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?

Yes, 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.

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)?

• 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.
• 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.

5. In Table 2: Time Performance, what's the time unit?

The time reported is in Seconds (s)

Dear reviewer grLC, we wanted to check if there are any remaining questions or concerns about our rebuttal. Since the author-reviewer discussion period ends tomorrow, we’d appreciate any final feedback you might have.

Thank you for your time and attention.

Dear Reviewer,

As the discussion deadline approaches, we are eager to address any remaining concerns you may have. We kindly request you to review our responses. If there is anything further we can do to facilitate the evaluation process, please let us know. We also respectfully ask that you reconsider our score based on the rebuttal provided.

Thank you for your attention.

Best regards,

Authors

I thank the authors for the detailed response. Overall, I think this is a solid paper and I would like to keep my rating of "accept".

Dear Authors,

Yes, that is appropriate.

Best,
AC