Learnable Stability-Aware Unstructured Grid Coarsening Using Graph Neural Networks for Accelerated Physics Simulations

Thank you for your response. I think my current score and confidence are appropriate.

Dear Reviewer,

Thank you very much for your mostly positive feedback. We look forward to providing comprehensive responses to your questions.

Q1. Are alternative remeshing strategies, such as Delaunay triangulation, considered, and if so, how do their performances compare?

Thank you for the question. We already use the Delaunay triangulation as a precondition for constructing the Voronoi tessellation in order to discretize the domain of the problem. There is a duality between a Voronoi tessellation and a Delaunay triangulation. We construct an updated Delaunay/Voronoi tessellation at every step of optimization. The finite volume method we employ for spatial discretization benefits from the Voronoi-based grid due to its flexibility in handling irregular geometries and ensuring local conservation properties.

Q2. Is your proposed method adaptable to non-fixed mesh topology systems, such as Lagrangian systems like the Deforming Plate model discussed in MeshGraphNets[1]?

Good question. As we observe, the mesh in the Deforming Plate example is a surface triangular mesh. The framework we propose is suitable for such meshes unless we have a mesh-based geometrical pooling algorithm: based on fuzzy clustering this algorithm should perform geometrical pooling but preserve the edges of the initial 3D body. Also, we should have a differentiable quasi-static simulator to pass gradients through. Moreover, if we substitute the mesh for another volumetric mesh, say 3D Voronoi volumes, then our method is also applicable. We currently work on 3D differentiable Voronoi tessellation.

Q3. Can you illustrate the meaning of colorbar in Figure 8 of Appendix A? Thank you for your comment! It is the value of the permeability. We will clarify it in the final version of the manuscript.

Dear Reviewer vcgH. We have incorporated new experiments and provided a detailed comparison of our algorithm with MeshGraphNets and other baselines. Additionally, we have thoroughly revised our paper, as outlined in our general replies. Our experiments reveal that surrogate-based approaches like MeshGraphNets, even when trained and inferred using the same data and timesteps as our method, exhibit lower inference accuracy. Furthermore, our method demonstrates superior data efficiency and our simulator works on coarse grid. We think we can integrate all remaining concerns due to camera-ready period so we kindly ask you to consider increasing your final score.

Q1. How does the GNN decides which nodes to keep in corsened mesh? Top-K?

In our experiments, we set which nodes to keep as a hyperparameter. The motivation for this problem comes from the field of proxy modeling and other engineering fields where model identification needs to be performed based on limited sensor measurements. For example, consider an oil reservoir with wells where measurements are available - we are only interested in simulating the pressure at these particular points. Thus, if we achieve speed-up of computation at these positions while preserving the physics, we consider the goal to be completed.

Indeed, if higher accuracy is required in certain regions of interest, additional sensor points can be incorporated into these areas to reduce interpolation errors. It can be implemented by the a priori analysis of the initial data: we may compute the gradient field and choose points that have higher gradients as sensor points. Also, we may use other algorithms based on graph spectral decomposition and other importance metrics e.g. centrality measures. This adaptive selection of sensor points ensures that regions requiring more precise solutions receive adequate attention without significantly increasing the computational gains.

Q2. Once coarsened nodes are obtained, how to get the new mesh? The current method only mentions triangular mesh, what is your proposal to adapt in these cases where the mesh is 3D or no longer simplex?

Thanks for pointing out this. Once the coarsened nodes are obtained, we build a new Delaunay mesh -> new differentiable Voronoi tessellation. This is a $nlogn$ operation where $n$ is not large since n is the number of the coarsened points ($n \ll N$). If we consider a 3D domain when the only change we should perform is to change the 2D differential Voronoi for a 3D setup (we are currently working on it). GNN and the Finite Volume simulator are agnostic to the 2D/3D embedding of the computational domain/graph. If the mesh is no longer simplex when it becomes quite hard to keep it regular after performing the pooling step, that is why we use Finite Volume solver as it allows a lot of flexibility. However, in general, once we get a differentiable coarsening and pooling for a regular grid our framework is perfectly suitable.

Q3. How many timesteps do you keep the gradients, 1, eg only calculating the gradients for stepping from n to n+1, or multi-step; if multi-step, how many?

We use a multi-step approach and keep gradients for the whole rollout (1000 timesteps).

Q4. You mentioned implicit method can improve results. But where are the results to support this claim?

Good question. We used the implicit solver to support the stability of our method because of the well-known stability properties of implicit solvers. We supported this claim by demonstrating the work of the implicit solver for the diffusion equation. The main advantage comes from the fact that for diffusion equation without Stability Loss explicit solver diverges. So, to broaden the scope of our method we implemented implicit solver as a building block for our framework. We provide the comparison of the solution for $\tau=0.01$ for the following scenario:

Setup

We set 100 time steps, with Gaussian-like initial condition for pressure with no source terms. And the comparison for our solvers:

Comparison of explicit and implicit JAX+Torch solver

By implementing the implicit solver we broaden the applicability of our method.

Q5. The simulation is required to be differentiable. This easily achieved in explicit scheme, but challenging to general PDEs solution schemes, eg, WENO for nonlinear advection terms. What is your proposal to these extensions?

Although there are several exceptions, explicit schemes are limited to application to non-stiff problems. So the question is basically how to design a differentiable solver to stiff problems, including nonlinear ones. Our vision is that differentiable implicit solvers, (backward Euler and more complex) can address this type of problem. It is one of the reasons why we have also considered differentiable implicit solver in our work.

Q6. In sec 3.3., you mentioned Sij is learnable by NN. What is the NN that outputs Sij? as Sij depends on the topology, ie the graph G, it should be output via some GNN strucutre? Otherwise it's not possible to extend to other mesh structure.

Thank you. As described in section 3.2 GRAPH COARSENING PIPELINE our NN architecture consists of a one graph convolution layer followed by a point-wise MLP with two fully connected layers, which predicts the cluster assignments Sij. We use softmax at the end. So, indeed, our architecture depends on the topology of computational graph. But we will improve our text for better readability.

Dear Reviewer,

Thank you for your detailed questions! We are glad to provide the answers below.

A very important point to note for all of these concerns is the different nature of [Shumilin et al.] baseline from usual baselines. In the [Shumilin et al.] is solved a very difficult problem of global optimization of mesh points, which is performed for every new grid. In the proposed by us approach the GNN outputs directly a solution of this optimization problem which is much more simple and general. That's why this baseline can be considered as the one we would like to be close with (so it's more about the best achievable metrics than about the metrics we can beat). So we consider a major strength of our approach that it gives results very close to the one by [Shumilin et al.].

C1. Insufficient Performance Advantages

Despite of the specific nature of baseline, we conducted a more thorough hyperparameter search. After that, we compared the effectiveness of our method with the two baselines we used in our paper. We considered the ”loop scenario”, this experiment was described in the paper in section 4.2. After selecting hyperparameters, more comparable results can be obtained with baseline of [Shumilin et.al.]. Please, see the detailed description of the experiment in our answer to W4 of Reviewer hygj.

C2. Limited Generalization Capabilities

It is important that our method does not need retraining from scratch for every new grid but the method from [Shumilin et.al.] requires retraining. When we mentioned improved generalization ability, we meant exactly this fact: our method has a trainable coarsening procedure that obtains a coarsened point cloud from the original one (and calculates corresponding features). The baseline method [Shumilin et.al.] does not have this procedure, it only optimizes the locations of points by running the overall optimization cycle from scratch for every new example. Our method is a great novel extension of this baseline with trainable coarsening, stability consideration and possibility of using implicit time-stepping inside the coarsening pipeline.

Additionally, we improved our trained model for performing grid coarsening at a reduction degree of 7.5% by adding more sensor nodes. Here the result of the training process with the result for pressure in zero sink.

Anonymous link: setup of the experiment with multiple sensor points

Next, we take the original permeability distribution field containing 312 points and find the coarsened point cloud using k-means + averaging approach, the resulting values serve as the initial permeability field for future comparisons. Here we introduce a new grid which contains 100 points and show our new results:

Anonymous link: new result

Here is the previous result from the paper:

Anonymous link: old result

We can see that only adding more sensor points improves the generalization ability of our model.

C3. Limited Computational Efficiency Improvement

Good question. A little misunderstanding occurred. According to the time comparison table in the paper, our method is slower than the baseline but of a similar order of magnitude. We justify this feature by the possibility of further reuse of the trained GNN. The baseline method should be run again if we want to apply it. We consider this as a big advantage of our method as we may exploit this feature when coarsening should be done repeatedly, for example in proxy modeling or inverse tasks.

Dear Reviewer WkfS, thank you for thorough review of our paper. With our explanations, we hope we have now addressed your concerns. If you have any additional concerns, please do let us know. We will be happy to address your comments. Otherwise, we kindly request you to support our work by raising the score. Additionally, we incorporated most of your suggestions about clarity of writing into our manuscript and continue incorporating remaining ones shortly.

Dear reviewer WkfS,

C3 not resolved; could you provide the simulation for original fine mesh that achieves similar accuracy? that can reveal the real efficiency gain Thank you for this insightful comment. For the simple loop scenario, we achieved similar accuracy on the coarsened grid (5%) as on the original fine grid (2x time speedup). For large scale scenarios we can achieve 10+x speedups. We appreciate the importance of providing analyzing computational gains, and we will ensure that this information is clearly highlighted in the camera-ready version of our manuscript.

Q1 not resolved; in many general cases, the "interesting point" will evolve over time.

Thank you for highlighting this important point. In our current framework, we focus on scenarios where the "interesting points" (e.g., sources, sinks, or key nodes) remain static over time. This is a common assumption in many physics-based simulations, where fixed key points such as wells in subsurface flow or boundary conditions in fluid dynamics are predefined and remain constant. However, we agree that in more general cases, "interesting points" could evolve dynamically over time. Extending our framework to handle such scenarios would involve integrating mechanisms to dynamically track and update these key points during the simulation. While this is beyond the current scope of our work, we recognize its importance for broader applications and plan to explore these extensions in future work. We will revise the manuscript in camera-ready version to acknowledge this limitation and outline potential directions for addressing it.

Q3 If you keep gradients for 1000 steps, will not that causes to anything like gradient diminishing or exploding? Can you clarify this? Thank you

Thank you for raising this important point. To prevent exploding gradients during long rollouts, we implemented gradient clipping, which ensures that gradients remain within a safe range. While gradient clipping does not directly address diminishing gradients, we have not observed issues with vanishing gradients in our experiments. We will clarify this in the manuscript in camera-ready version to address your concern.

Additionally, please see our general replies about comparison with MeshGraphNets and other baselines and all modifications we have made in our paper already. I think we can integrate all remaining concerns due to camera-ready period so we kindly ask you to consider increasing your final score.

Q1. "The largest question I have is that, since a graph neural network is already employed, why not further the effort and directly generate the solution with the graph neural network (by increasing the current one layer network to a slightly larger one)?"

We view the use of GNNs for coarsening and as a surrogate solvers as complementary rather than mutually exclusive. Each approach has its advantages:

• GNN for Coarsening + Differentiable Numerical Solver. Pros: Maintains high physical accuracy and fidelity, generalizable to various scenarios, ensures solutions adhere to physical laws. Cons: May be slower than surrogate models due to the computational cost of the numerical solver.
• GNN as Surrogate Solver. Pros: Faster inference times, potential for real-time simulations. Cons: May sacrifice accuracy, less generalizable, requires extensive training data for different scenarios, leads to exploding errors when run on a long trajectory.

By focusing on coarsening, we aim to reduce computational costs while retaining the reliability of numerical solvers. This approach is particularly beneficial when high accuracy is required at specific critical points or under varying physical conditions.

While directly generating the solution with an expanded GNN-based surrogate model is an intriguing idea, our current approach prioritizes maintaining accuracy respecting the physical laws and generalizability through the use of numerical solvers after coarsening. We believe that this strategy effectively balances computational efficiency with the need for accurate, physically consistent solutions. We appreciate your thoughtful suggestion and are excited about the possibility of exploring this avenue in future research.

Q2. "Table 1 & 2, it is reported that the method is actually slower and more memory intensive than the baseline, while results are comparable or only slightly better???"

Thank you for noting this. Indeed, the method is slower and memory intensive than the baseline. However, this cost is needed to make the method adaptive and generalizable without re-training, whereas the baseline needs re-training for each task. Additionally, please see the answer to the reviewer hygj where we obtained better hyperparameters.

Dear Reviewer, Thank you for your feedback. We appreciate you pointing out the strengths of our work. We address all of your questions in the following.

W1. "Whilst the solution from the coarsened grid are accurate at the node points, the super-coarsened grid inevitably lead to another source error since quite a great proportion of applications would require the field solution rather than discrete points, which would necessitate interpolations from the solution points..."

We appreciate your insightful comments regarding the potential interpolation errors when reconstructing the full pressure field from the coarsened grid. We want to highlight that in our work, we focus on control (sensor) points and do not assume access to the full field solution. Our objective is to develop a coarsened grid (to reduce the computational cost) while minimizing the error on the sensor points where measurements are available. The motivation for this problem comes from the field of proxy modeling and other engineering fields where model identification needs to be performed based on limited sensor measurements. For example, consider an oil reservoir with wells where measurements are available - we are only interested in simulating the pressure at these particular points. Thus, if we achieve speed-up of computation at these positions while preserving the physics, we consider the goal to be completed.

Indeed, if higher accuracy is required in certain regions of interest, additional sensor points can be incorporated into these areas to reduce interpolation errors. It can be implemented by the a priori analysis of the initial data: we may compute the gradient field and choose points that have higher gradients as sensor points. This adaptive selection of sensor points ensures that regions requiring more precise solutions receive adequate attention without significantly increasing the computational gains.

Multiple sensor points setup

Multiple sensor points setup

To demonstrate the effectiveness of our approach in predictive tasks at specific measurement points, we conducted experiments modeling predictions for future time steps not used during the training stage. In the "loop" scenario with 8 sinks (representing critical points), we compared our method against the baseline approach by [Shumilin et al.]. Importantly, the training loop is performed only on 1000 timesteps for both algorithms.

Prediction experiment 0

Prediction experiment 1

Prediction experiment 2

Prediction experiment 3

Results show that our algorithm beats the competing algorithm [Shumilin et.al.] for sinks 0, 1, 2, 4, and 6. Hence, it demonstrates competitive or better performance and, additionally, does not only optimize points but represents a trainable coarsening approach, providing a novel extension to the algorithm by [Shumilin et.al.]. Possibly, a better choice of hyperparameters could lead to even better results.

Altogether, our algorithm coarsens the grid not using the solution over all the grid but the solution only at several critical points.

To clarify these features of our method, we will include the additional experiment above and new text in the original manuscript.

Thank you for your responses, which clarify a series of issues, but unless other reviewers deem the work good enough, I still cannot suggest the work to be published in ICLR. I have raised the score to 5, but I have to check the comments from other reviewers before deciding whether to increase it further to 6.

Dear reviewer wVwb, please see our general replies where we compared our algorithm with other baselines and have significantly improved our manuscript. We kindly ask you to increase your final score considering the work we have done.

Dear reviewer, thank you for your thoughtful feedback. We appreciate your efforts. We provide answers below:

W1. “They claimed the differentiable implicit method is their novelty; however, PyTorch supports differentiable linear solvers”

Answer for W1. We have not claimed an implicit method as our novelty, we have claimed an implicit differentiable Finite Volume solver to be used in our coarsening framework. However, we agree that our current formulation contribution sounds misleading and we will correct the formulation in our manuscript as follows: We develop the differentiable implicit numerical finite volume solver that leverages differentiable Voronoi tesselation.

While PyTorch’ torch.linalg.solve is a powerful tool for dense linear systems, it lacks native support for efficient sparse matrix operations. In contrast, JAX, particularly with its spsolve method from jax.experimental branch is optimized for handling large-scale sparse linear systems. This is crucial for real-world tasks where the underlying matrices are inherently sparse due to the nature of finite volume discretization.

W2. “In addition, the community knows that there are implicit GNNs that solve implicit problems”

Answer for W2. Yes, implicit GNN is a known approach. However, implicit GNNs primarily focus on enhancing message passing mechanisms to capture long-range dependencies within graphs without being constrained by a fixed number of layers. The term “implicit” in this paper means using implicit updating scheme for message passing instead of a classical explicit updating message passing mechanism. This is completely different from our “implicit” meaning in the context of implicit scheme for solving linear systems.

W3. “but there are actually a lot of works investigating mesh coarsening”

Answer for W3. We would like to point out that none of the works similar to ours have the implicit differentiable finite volume solver without any approximations (non-surrogate):

[1] https://arxiv.org/pdf/2405.04466 (A fully differentiable GNN-based PDE Solver: With Applications to Poisson and Navier-Stokes Equations). Authors developed a differentiable finite volume approach using FVM loss for GNNs. They use the Encoder-Decoder approach with GNN blocks inside. However, they use PDE loss for training. Therefore, the solver they use can be considered as a surrogate.

[2] https://arxiv.org/pdf/2010.03409 (LEARNING MESH-BASED SIMULATION WITH GRAPH NETWORKS)

MeshGraphNet applies remeshing which does not imply coarsening, in contrast, it may significantly increase the number of cells and load the computations. In fact, the remesher in MeshGraphNets coarsens the mesh by removing edges when removing an edge does not create any invalid edges; the criterion for invalidness is $u_{ij}^T S_{ij} u_{ij} \leq 1$, where $u_{ij}$ is the vector representation of the edge $ij$, $S_{ij}$ is the sizing field tensor which is predicted. This is a purely geometrical criterion for coarsening, in contrast, our method uses simulation data misfit to cluster the points and further aggregate them. Moreover, MeshGraphNet remesher is only applied to the 2D manifolds embedded into 3D like Flag or Sphere scenarios as it's based on the local curvature.

[3] https://arxiv.org/abs/1904.09019 (Graph Element Networks: adaptive, structured computation and memory) GENs approximate the solution of PDEs using GNN as a surrogate model, trained on high-resolution solutions. They do not integrate the actual physics solvers into the training loop. GENs adjust node positions or densities to better capture the solution space but do not perform explicit graph coarsening through node aggregation or clustering. Additionally, Graph element networks can solve only the stationary partial derivative equation in their implementation or predict displacement. The authors have not included in the code the possibility of taking into account time dependence.

W4. “The experimental evaluation is weak.” Answer for W4. We compared the effectiveness of the three methods at different levels of grid reduction. Consider the ”loop scenario”, this experiment was described in the article in section 4.2. After selecting hyperparameters, comparable results can be obtained. The following hyperparameters were used for the reduction degree of 10%:

• Time Step: 0.0001
• Number of Epochs: 300
• Learning Rate: 0.015
• Sigmoid Weight: 10
• Physics Loss Weight: 15
• For reduction degree 7.5%:
• Time Step: 0.0001
• Number of Epochs: 300
• Learning Rate: 0.01
• Sigmoid Weight: 10
• Physics Loss Weight: 20
• For a reduction degree of 5%:
• Time Step: 0.0001
• Number of Epochs: 300
• Learning Rate: 0.01
• Sigmoid Weight: 10
• Physics Loss Weight: 15

Experiment results 1: anonymous link

Next, consider the case of the Sinusoidal variable permeability scenario. Let's demonstrate the experiment described in section 4.3 in which we compare RMSE for three algorithms in two different sinks. It can also be seen here that after selecting hyperparameters, we got a comparable result. The following hyperparameters were used for reduction degree 10%:

• Time Step: 0.0001
• Number of Epochs: 300
• Learning Rate: 0.015
• Sigmoid Weight: 10
• Physics Loss Weight: 20
• For reduction degree 7.5%:
• Time Step: 0.0001
• Number of Epochs: 300
• Learning Rate: 0.01
• Sigmoid Weight: 10
• Physics Loss Weight: 20
• For reduction degree 5%:
• Time Step: 0.0001
• Number of Epochs: 300
• Learning Rate: 0.01
• Sigmoid Weight: 10
• Physics Loss Weight: 25

Experiment results 2: anonymous link

Experiment results 3: anonymous link

Additionally, we added some preliminary experiments with adaptive mesh refinement. See our answer to the question below.

W5. “The effectiveness of the proposed method is limited. The experimental results suggest that the method has almost the same accuracy and longer computation time for training. If the reviewer did not miss, there is no evaluation of the prediction time, which is one of the most essential metrics for the mesh coarsening task.”

A very important point to note is the different nature of [Shumilin et al., 2024] baseline from usual baselines. In the [Shumilin et al., 2024] is solved a very difficult problem of global optimization of mesh points, which is performed for every new grid. In the proposed by us approach the GNN outputs directly a solution of this optimization problem which is much more simple and general. That's why this baseline can be considered as the one we would like to be close with (so it's more about the best achievable metrics than about the metrics we can beat). So we consider a major strength of our approach that it gives results very close to the one by [Shumilin et al., 2024].

As for the prediction time, let us clarify what the prediction time means in this context. For the baseline from [Shumilin et al., 2024], this time is the training time from scratch, since this algorithm cannot run without running the optimization from scratch. However, our algorithm can be trained once and predicted on new examples without re-training. Obviously, its prediction time does not include training from scratch and this distinguishes our algorithm from the baseline [Shumilin et al., 2024].

We also use a comparison with the coarsening algorithm from METIS. METIS is a software tool designed for partitioning graphs and meshes. METIS aims to generate blocks that are uniform in size and shape, minimizing the number of connections between different partitions. If we let fine-scale transmissibility measure connection strengths between cells for the edge-cut minimization algorithm, the software tries to construct a block without crossing large permeability contrasts. We made a comparison using the examples described in our paper. First result for Sinusoidal variable permeability scenario:

Experiment results 4: anonymous link

The second result for the “loop” scenario:

Experiment results 5: anonymous link

The result obtained by our algorithm is much more accurate than the result of METIS .

Dear Reviewer hygj. Thank you for your comment! We have added new experiments and comparison of our algorithm with MeshGraphNets and FluxGNN. Also we have thoroughly updated our paper. See our general replies. Experiment shows that such surrogate-based approaches as MeshGraphNets, even if we use same data, same timestep for training and same timestep for inference as our method, demonstrate worse accuracy in inference even on fine grid. But our method is more data-efficient and our simulator works on coarse grid (obviosly, because our algorithm is coarsening algorithm). Additionally, we checked FluxGNN and their code is irreproducible (see also our general reply).

We thank you for thoughtful comment and agree that implicit GNNs and implicit numerical methods share a conceptual foundation, as both involve solving algebraic equations iteratively. While the goals and areas of application are somewhat different, it would indeed be interesting to explore the possibility of designing a surrogate solver based on implicit GNNs in future work. However, in this study, we focus on integrating an exact numerical solver using implicit schemes (e.g., backward Euler) to ensure precise and stable solutions for physical simulations without relying on approximations or surrogate models. Also we developed and used fully diffferentiable Voronoi tessellation. We will add a clarification in the camera-ready to highlight this distinction and the unique focus of our work.

Additionally, we will add analysis of speed-accuracy tradeoff in some scenarios while camera-ready period.