{'fraction_non_alphanumeric': 0.030348000296801958, 'fraction_numerical': 0.010388068561252505, 'mean_word_length': 4.751013441433753, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
M2N: Mesh Movement Networks for PDE Solvers

Mainstream numerical Partial Differential Equation (PDE) solvers require discretizing the physical domain using a mesh. Mesh movement methods aim to improve the accuracy of the numerical solution by increasing mesh resolution where the solution is not well-resolved, whilst reducing unnecessary resolution elsewhere. However, mesh movement methods, such as the Monge-Ampere method, require the solution of auxiliary equations, which can be extremely expensive especially when the mesh is adapted frequently. In this paper, we propose to our best knowledge the first learning-based end-to-end mesh movement framework for PDE solvers. Key requirements of learning-based mesh movement methods are alleviating mesh tangling, boundary consistency, and generalization to mesh with different resolutions. To achieve these goals, we introduce the neural spline model and the graph attention network (GAT) into our models respectively. While the Neural-Spline based model provides more flexibility for large deformation, the GAT based model can handle domains with more complicated shapes and is better at performing delicate local deformation. We validate our methods on stationary and time-dependent, linear and non-linear equations, as well as regularly and irregularly shaped domains. Compared to the traditional Monge-Ampere method, our approach can greatly accelerate the mesh adaptation process, whilst achieving comparable numerical error reduction.


Introduction
Partial Differential Equations (PDEs) are widely used to model natural phenomena, ranging from aerospace engineering, atmosphere/ocean dynamics, bio-engineering to computer graphics. Efficient and accurate numerical solutions for complex PDEs are essential but challenging problems in all scientific and engineering disciplines.
Solving PDEs using numerical methods such as the finite element method (FEM) requires discretizing the problem spatially and temporally. The spatial discretization is usually represented using a mesh. The quality of a mesh affects the accuracy of the numerical solution. However, it is often prohibitively expensive to solve the problem on a very high resolution mesh. Mesh adaptation is an advanced discretization method to tackle this problem. It increases the mesh resolution where the solution requires higher numerical accuracy, while decreasing the mesh resolution where unnecessary. Mesh adaptation methods can be generally divided into two categories: h-adaptation and r-adaptation. In h-adaptation, new mesh nodes are dynamically added to the regions where fine resolution is required. In r-adaptation (or mesh movement), mesh nodes are only relocated or moved without changing the mesh topology. Compared to h-adaptation, r-adaptation has several attractive features. First, no extra mesh points are generated, which keeps the dimension of the linear system representing the discretized PDE unchanged. In addition, fixed mesh connectivity can also make the structure of the stiffness matrix unchanged, which enables matrix pre-factorization that can accelerate the solution of the large linear systems encountered in FEM ]. However, a key challenge for mesh movement methods is mesh tangling, whereby one or more elements become inverted. Mesh movement methods based on optimal transport theory [Weller et al., 2016, McRae et al., 2018 can effectively prevent mesh tangling issues (see [Clare et al., 2022] for a discussion on this), but require solving a Monge-Ampère equation at each adaptation step, which is extremely time-inefficient.
AI for the solution of PDEs has been an emerging topic in recent years, and shows great potential in solving problems where traditional numerical PDE solvers struggle (e.g. high dimensional problems [Han et al., 2018, Sheng and), or in accelerating the solution process by learning a neural operator from the parameterized description of a PDE problem to its corresponding solution [Li et al., 2020, Lu et al., 2019. However, these methods still encounter fatal bottlenecks, such as the precision guarantee of the solution, data efficiency, generalization capability, and so on. These are fundamental limits of deep learning, but essential for the scientific computing scenarios.
A possible alternative is to perform learning based mesh adaptation, by which the traditional numerical PDE solvers can achieve better performance, while the time consumption of mesh adaptation is greatly reduced. There have been some works in this direction [Zhang et al., 2020, Fidkowski and Chen, 2021, Tingfan et al., 2021, Pfaff et al., 2020. However, most previous methods target mesh generation [Zhang et al., 2020] or mesh refinement , Fidkowski and Chen, 2021, instead of topology-invariant mesh movement as considered here. Moreover, the previous methods are not end-to-end, which means the neural networks are used to predict certain metrics, such as the local mesh density [Zhang et al., 2020 or the metric tensor [Fidkowski and Chen, 2021, Tingfan et al., 2021, Pfaff et al., 2020, which then have to be fed into a traditional mesher/remesher to obtain the mesh. Therefore the overall performance is bounded by the mesher/remesher, which is computational geometry based and hence not optimized for solving PDE problems. On the contrary, in our method, the adapted mesh is directly output by the neural network.
In this work, we propose to the best of our knowledge the first learning-based end-to-end mesh movement framework for PDE solvers. Taking the source term, the input field, and/or the PDE parameters as input features, the model deforms an initial mesh to the adapted mesh by mesh movement. In usage, the model can be applied to a class of PDEs without retraining. We design a Neural-Spline based model for mesh deformation. It is an invertible neural network and hence can avoid mesh tangling. Moreover, its mechanism naturally guarantees a hypercubic boundary can be maintained through the learnable mapping. We also design a graph attention network (GAT) based model for mesh deformation. The graph neural network can naturally describe domains with irregular shapes and embed the relevant information. We utilize the attention mechanism of the GAT model to guarantee each mesh node stays within its neighborhood.
Our main contributions can be listed as below: 1. We propose a learning-based end-to-end mesh movement framework for PDE solvers, which to the best of our knowledge is the first of its kind. Without interfering with the PDE solver, the models can achieve numerical error reduction similar to the traditional Monge-Ampère method, while the mesh generation process is accelerated by two to three orders of magnitude.

2.
A Neural-Spline based model and a GAT based model are proposed for mesh deformation. Besides generalization to different PDE parameters, source terms, input solution field, etc., the models are specifically designed to guarantee boundary consistency, alleviate mesh tangling, and generalize to different mesh densities.

Related Work
Mesh movement method. Mesh movement methods include velocity-based and location-based methods. In the work of [Anderson and Rai, 1983, Gnoffo, 1982, Farhat et al., 1998], the mesh is moved according to attraction and repulsion pseudo-forces between nodes motivated by a spring model. The moving mesh finite element method [Baines et al., 2005] computes the solution and the mesh simultaneously by minimizing the residual of the PDEs written in a finite element form. As for location-based method, the moving mesh PDE (MMPDE) method [Huang and Russell, 2011] moves the mesh through the gradient flow equation of an adaptation functional. In recent years, there has been a growing interest in optimally-transported r-adapted meshes [Budd and Williams., 2009, Weller et al., 2016, McRae et al., 2018, Clare et al., 2022.
AI for PDE. To solve a PDE problem, neural networks can be used to represent the function to solve, and trained either with the residual loss of the PDE or using the variational principle [Raissi et al., 2019, Yu et al., 2017. Neural operators are proposed to learn an operator from the problem function to the solution function [Li et al., 2020, Lu et al., 2019. There also exist mesh-based PDE solvers with deep learning. In [Pfaff et al., 2020], a graph neural network with additional world edges are applied to predict dynamical systems, shown to be effective with a wide range of physical systems. In [Belbute-Peres et al., 2020], a differentiable PDE solver is embedded in a neural network to help predict accurate solutions and also backpropagate the loss so that the input coarse mesh can be optimized.
AI for Meshing. AI methods have also been proposed for mesh generation, adaptation, and so on. MeshingNet [Zhang et al., 2020] uses a neural network to learn the required local mesh density, which can then be provided to a Delaunay triangulation based mesh generator to generate high quality meshes. The optimal local mesh density is also learned in  for mesh refinement. On the other hand, the mesh refinement process is formulated as a reinforcement learning problem in  to minimize the PDE solution error under given refinement budgets. The flow field is predicted by machine learning models to calculate the metric tensor so that the mesh can be optimized accordingly [Tingfan et al., 2021]. The authors in [Fidkowski and Chen, 2021] focused on optimal anisotropic meshes by predicting the desired element aspect ratio. In [Pfaff et al., 2020], the sizing field is predicted by a neural network for adaptive remeshing along with the system dynamics.

Problem Statement
Meshing is the procedure to spatially discretize a PDE, which is necessary for most numerical methods to solve PDEs. Among these methods, the Finite Element Method (FEM) is one that has been widely used in various engineering fields.
A high-quality mesh can significantly improve the accuracy-efficiency trade-off of PDE numerical methods, therefore mesh adaptation techniques have been well studied, which in major include two categories, mesh refinement and mesh movement. In this work, we apply our mesh movement method with FEM.
Consider a discretized PDE P defined on the domain Ω in D dimensions. For mesh adaptation, we can define two domains: Ω i for the initial mesh T i , and Ω d for the adapted mesh T d . T d should be adapted in response to the structure of the physical system under consideration. We denote mesh node positions in T i by ξ, and in T d by x. The adapted mesh T d is the image of T i under the mapping x = F(X , ξ|θ), where X is the input state of the PDE and θ represents the model parameters. The solution on T i is denoted as u i = P(X (ξ)) serving as a baseline and the solution u d = P(X (x)) on T d is the optimized solution enhanced by mesh adaptation. The computational cost of traditional methods for mesh adaptation is often expensive. For some methods, the cost of mesh adaptation is comparable to or even more expensive than that of solving the underlying PDE problem, which is generally unacceptable. Therefore, our goal is to model the mapping F(·|θ) using learning based methods such that the mesh adaptation process can be greatly accelerated.

Framework Overview
The proposed framework is demonstrated in Figure 1(a). For a class of PDEs, each of its specific sample P can be parameterized by PDE related information, such as source terms, boundary conditions, PDE parameters, solution fields, etc. We use X to denote such information and take it as the input of the neural network, so that the trained model can generalize to a class of PDEs. The mesh movement network F(·|θ) deforms the initial mesh T i and outputs an adapted mesh T d , given the input state X . T d is then fed into the PDE solver to help improve the accuracy of the solution.
Under the M2N framework, we consider two implementations of the mesh movement network: a Neural-Spline based network and a GAT based network, which are denoted as M2N-Spline and M2N-GAT, respectively.

PDE Solver
Initial Mesh

Input State
Adapted Mesh

Output State
Mesh Movement Network

GNN Block GAT Deformer
Global Feature Extractor

Graph Constructor
Physical Parameters

Sample
Conv GAP

Neural Spline
Global Feature Extractor

Neural-Spline based Network
Our Neural-Spline based network is mainly composed of two parts. One is a global feature extractor to extract features from the input state X , the other is a neural spline deformer to output the adapted mesh T d .
Global Feature Extractor We uniformly sample the input state X in the domain Ω i . If the domain is with an irregular boundary, sampling is performed inside its minimum bounding box. The values of the sampling points outside the domain boundary are set to zero. The sampled data are assembled as a state tensor M. To keep the extracted feature invariant to the magnitude, the state tensor M is normalized by its maximum absolute value. After normalization, M is sent into the convolutional layers for feature extraction, whose output is further fed into a global average pooling layer to obtain a global embedding E, which is invariant to the mesh resolution. The embedding E is then concatenated with physical parameters such as the mesh density and the PDE coefficient to get I, which is fed into the deformer.
Neural-Spline based Deformer Normalizing flow models [Kobyzev et al., 2020] are proposed to learn invertible mappings. Neural spline [Durkan et al., 2019], as a specific type of normalizing flows, transforms each dimension of the input with a differentiable monotone rational-quadratic spline function S(·). More specifically speaking, for the d-th dimension of the node coordinates ξ of the initial mesh, the anchor points of the spline function are parameterized by the input features I and the coordinates of the other dimensions ξ −d , to get the output S d (ξ d |I, ξ −d ). In a neural spline block, each dimension takes turns to be transformed, so that the node coordinates ξ of the initial mesh will be mapped to node coordinates x to assemble the adapted mesh T d .
Since the neural spline model is guaranteed to be an invertible mapping, mesh tangling can be alleviated. A specialty of the neural spline model is that, the end points of the spline function are fixed, therefore the intervals of the input and the output can be kept unchanged. In our case, this property is utilized to preserve the mesh with a hypercubic boundary (e.g. a rectangle in the 2-D case). Moreover, since the neural spline learns a continuous mapping, it can naturally be generalized to different mesh densities.

GAT based Network
Although the Neural-Spline based network is well-designed for meshes of hypercubic domains, it is difficult to extend to domains with more general boundaries. Therefore, we propose a more flexible model based on the Graph Attention Network (GAT) [Veličković et al., 2017]. It consists of a feature extractor and a GAT based mesh deformer. Figure 2: The implementation of the GAT Deformer. ξ n and H n denote the mesh node positions and extracted features for the nth GAT block, while H n represent the corresponding outputs.

Feature Extractor
The feature extractor of the GAT-based network consists of two branches: one is the Global Feature Extractor which is similar to the Neural-Spline based network, the other is a Graph Neural Network (GNN) branch for local feature extracting.
For the GNN branch, a graph G = (V, E) is constructed based on the initial mesh T i . Mesh nodes and mesh edges are transformed to graph nodes V and bidirectional graph edges E. Then, we sample the values of the input state X at the position of each corresponding graph node v on the initial mesh. To capture the interaction between mesh nodes, the sampled values and mesh nodes geometry information are further processed by a GNN Block with message passing mechanism which can propagate the sampled physical quantities across the graph. To improve the generalization ability on mesh with different densities, the displacement vector u ij = u i − u j and its norm u ij are selected to describe the mesh density, where u i denotes the position of mesh node i. They are encoded into mesh edges e ij ∈ E as the relative edge features. These edge features are then updated by a Multi-Layer Perceptron (MLP) f with three layers and the features of each mesh node v i are updated by summing up the surrounding edge features: The features of each mesh node are concatenated with the extracted global feature h g and physical parameters p (such as the viscosity coefficient in Burgers equation). The concatenated feature H together with the mesh nodes position ξ are sent into the GAT deformer for mesh movement. Figure 2, the deformer consists of N connected GAT blocks with SELU as the activation function. Each block takes two streams of inputs: mesh vertex positions ξ and extracted features H. To improve the learning efficiency, the intermediate output mesh nodes positions ξ n are appended to the corresponding output features H n as the part of input features H n+1 for next block. The GAT based deformer can naturally handle the problem of mesh tangling: during the deformation process, the movement of each mesh node is confined inside the convex hull composed of its 1-ring neighbors, which can effectively alleviate mesh tangling. To keep the boundary consistent before and after the deformation, the nodes on the boundary are restricted to move along the domain edge.

Experiment
We evaluate our proposed models against two baseline models using error reduction ratio, mesh generation speed, tangling avoidance, as well as generalization capability to different PDE problems and mesh densities.

Experiment Setup
Comparison Baseline. As mentioned in Section 1, there is no similar previous work to compare results against. Therefore, we compare our proposed models with two baselines that can be interpreted as ablated versions of our proposed models. The model MLP-Deform-Clip replaces the neural-spline block with MLP. The model GAT-Deform-Clip replaces the GAT-deformer with the ordinary GAT block. Instead of learning the positions of mesh nodes, we set the learning target as the mesh nodes displacement for the baseline models because it gives better performance. In order to enforce that baseline models can also preserve boundary consistency, the nodes moved out of the boundary will be pulled back into the domain, and the displacement component perpendicular to the boundary is clipped for the nodes which are supposed to stay on the boundary, as shown in Figure 3. However, there is no trivial way to alleviate mesh tangling. Monitor Function. In all experiments, we use the Monge-Ampère (MA) method as described in [McRae et al., 2018] to generate the target deformed mesh. For the MA method, the mesh is equidistributed with respect to a user-defined monitor function. The monitor function provides a concept of 'mesh density' and controls the sense in which errors are equidistributed by the mesh adaptation algorithm. Therefore, the choice of monitor function greatly impacts the geometry of the deformed mesh and should be chosen with care. In this paper, we define the monitor function as where u is the numerical solution on the current mesh, u exact is the numerical solution on a uniformly refined mesh or analytical solution, H(u) F represents the Frobenius norm of the Hessian matrix of the solution u.
Implementation Details. For all experiments, meshes are evaluated using PDE solutions obtained from the opensource FEM solver Firedrake [Rathgeber et al., 2016]. The neural network models are implemented with PyTorch [Paszke et al., 2019], and the graph neural networks are implemented with PyTorch Geometric (PyG) [Fey and Lenssen, 2019]. We use the 1 loss to train the models, i.e., L(ξ,x; θ) = F θ (ξ) −x 1 , wherex denotes the positions of target mesh nodes. The models are trained with the Adam optimizer [Kingma and Ba, 2014], using a mini-batch size of 6 meshes. Each model is trained three times with different random seeds.
Evaluation Metric. We evaluate the performance of different models by calculating the average error reduction ratio of the deformed mesh generated by different models. The quality of the mesh is evaluated by calculating the percentage reduction of the discretization error. The ground truth is provided either by the numerical solutions on meshes with very high resolution or the analytical solutions by proper dataset generation. We calculate the percentage of error reduction with the formula (e initial − e adapted )/e initial , where e initial is the error norm of the solution solved on initial uniform mesh compared with the ground truth. Similarly, e adapted is the error norm of the solution obtained from adapted deformed mesh.

Poisson's Equation
Poisson's equation is a second order, linear, stationary PDE, which is widely used in electrostatics and thermodynamics, amongst other fields. We consider the 2-D case with a Dirichlet boundary condition: −∇ · ∇u(x, y) = f (x, y), (x, y) ∈ Ω, u(x, y) = u 0 (x, y), (x, y) ∈ ∂Ω. (2) For both the square and heptagonal domain experiments, we generate analytical u samples from a mixed Gaussian distribution, which are fed into the Poisson's equation to obtain the corresponding source terms f and boundary conditions u 0 as the problem samples, whereas the u functions serve as the ground truth.
Square Domain In this experiment, we train the models on cases with mesh resolution of 15x15 and 20x20, each with 275 samples. To evaluate the models and how well they generalize to different mesh resolutions, we test on cases with mesh resolution from 12x12 to 23x23, each with 125 samples. Moreover, we deliberately set the optimal mesh movement to be drastic, in order to test how well different models can handle mesh tangling.  The quantitative results are summarized in Table 1. The proposed models, M2N-Spline and M2N-GAT, can achieve similar error reduction compared to the traditional MA method, while the mesh generation speed is two to three orders of magnitude faster. In addition, although the proposed models perform only marginally better than the two baseline models in error reduction and time, they are proven very effective to keep the mesh untangled. In comparison, the baseline models suffer from mesh tangling.
An example is given in Figure 4, where the mesh density in the upper left corner needs to be very high (shown in Figure  4(a)). It is demonstrated that, M2N-Spline is more flexible for the cases where the overall mesh deformation is required to be large. On the other hand, for M2N-GAT, because of the constrained movement in each layer to alleviate mesh tangling and the finite layer numbers, it does not learn as well in such scenarios.

Irregular Heptagonal Domain
To evaluate the performance of different models for more general domain shapes, we conduct an experiment using Poisson's equation in an irregular heptagonal domain. The models are trained at mesh densities of 13, 16, 19, and 22, each with 320 samples, and tested on mesh densities from 12 to 23, each with 80 samples, to evaluate the performance and generalization capability of the models.
The results are summarized in Table 2. It can be seen that all deep learning models can perform mesh movement around two to three orders of magnitude faster than the traditional MA method. Two GNN-based models perform better than the other two models, with the error reduction ratio comparable to the MA method. This is because the GNN-based models can naturally embed the information of the entire irregular domain into the network, and there are extra local feature extractors in the model. On the contrary, both the M2N-Spline and the MLP-Deform-Clip models are point-to-point mappings with only the global feature extractor, hence lack such capabilities. The results of all methods for an example problem are shown in Figure 5, from which we can see that the GNN-based models can better capture  the delicate local structure where mesh resolution needs to be increased. Although mesh tangling still will not occur for M2N-GAT, it happens for M2N-Spline, because for the M2N-Spline model, its guarantee only works for hypercubic boundaries (rectangle in the 2-D case).

Burgers' Equation
The viscous Burgers' equation is a non-linear, time-dependent PDE describing advection and diffusion processes in fluids. Here we consider the 2-D Burgers' equation, where constant scalar ν represents viscosity and u is the velocity vector field obeying this PDE. We define the problem on the unit square domain Ω = [0, 1] 2 .
In this experiment, we still train the models on cases with mesh resolution of 15x15 and 20x20, each with 9 trajectories and 60 timesteps per trajectory, with different viscosity coefficients. To evaluate the models and how well they can generalize to different mesh resolutions, we test on cases with mesh resolution from 11x11 to 24x24, each with 8 trajectories and 60 timesteps per trajectory.
The results are summarized in Table 3. It can be found that all learning-based methods are three to four orders of magnitude faster than the traditional MA method. The acceleration is about one order of magnitude larger than the  Poisson's equation experiments, because the MA method runs even slower for a nonlinear PDE. Mesh tangling never occurs for M2N-Spline and M2N-GAT in this experiment. The GNN-based models perform better than the other two models, because they contain extra local feature extractors, and local information can be propagated better through edges, which is important for the Burgers' equation problem. Some generated mesh examples at different timesteps of four different trajectories are shown in Figure 6, where it can be seen that M2N-GAT is better at performing delicate local deformation compared to the M2N-Spline model.

Conclusion
In this paper we have proposed the Mesh Movement Network (M2N), which to the best of our knowledge is the first learning-based end-to-end mesh movement method for PDE solvers. Traditional mesh movement methods can improve the accuracy of numerical PDE solutions without modifying the topology of the mesh, at the expenses of solving an auxiliary PDE, which is often very computationally expensive and sometimes makes the approach infeasible. With the power of deep learning, M2N generates adapted meshes for different PDE problems of the same type, with the solution precision comparable to ground truth but at a much faster speed. To achieve this robustly, we have designed a Neural-Spline based and a GAT based mesh deformer, to guarantee the output adapted mesh retains boundary consistency, alleviates mesh tangling, and generalizes to different mesh densities. The results are validated on the static linear Poisson's equation with regular and irregular domains, and the time-dependent nonlinear Burgers' equation.