PDENNEval: A Comprehensive Evaluation of Neural Network Methods for Solving PDEs

1 Introduction

Owing to their exceptional approximation capabilities [Hornik, 1991; Barron, 1994] and multi-level feature extraction abilities [LeCun et al., 2015], neural networks (NN) based methods have been widely applied and have achieved great success in various disciplines, including computer vision [Krizhevsky et al., 2012; He et al., 2016; Dosovitskiy et al., 2020] and natural language processing ([Vaswani et al., 2017; Devlin et al., 2018; Brown et al., 2020]). Recently, NN methods have excelled in solving partial differential equations (PDEs), yielding significant breakthroughs in this field, e.g. [E et al., 2017; E and Yu, 2018; Raissi et al., 2019; Li et al., 2020; Lu et al., 2021a]. Compared to traditional numerical PDEs methods like finite element methods, NN methods exhibit unique advantages. One notable advantage is that they are mesh-free, making them immune to the curse of dimensionality. In addition, NN methods can learn patterns and features from data which might be governed by multiple types of equations[Raissi and Karniadakis, 2018]. This capability enables a single solver to handle various types of equations, offering a more versatile approach for solving PDEs.

∗An extendede version of this paper(with the Appendix included) can be find in https://github.com/Sysuzqs/PDENNEval.

†Corresponding author.

The widespread development of NN methods for solving PDEs generates a pressing need for the assessment and comparative analysis of these techniques. In 2021, Lu and his colleagues compare the Physics-Informed Neural Networks (PINN) to a standard finite element method across various problems, including Poisson equation, Burgers' equation, Volterra integral differential equation, Lorenz system and Diffusion-Reaction systems in [Lu , 2021b ]. In 2022, three benchmarks of NN methods for PDEs are published. First, Lu and his colleagues[Lu et al., 2022] present the relative performance of Deep Operator Network (DeepONet) and Fourier Neural Operator (FNO) by using 16 different benchmarks covering Burgers' equation, Darcy problem, Advection equation, Compressible Euler equations, electroconvection problem and Navier-Stokes equation, etc. Secondly, a side-by-side analysis, called PDEArena, is presented in [Gupta and Brandstetter, 2022] to show the performance of NN methods FNO, ResNet, and U-Net on Shallow-Water equations and Navier-Stokes equations. Thirdly, a versatile benchmark suite[Takamoto et al., 2022b], called PDEBench, is proposed to compare NN methods FNO, U-Net, PINN, using datasets derived from 11 hydromechanical field PDEs, including Advection, Navier-Stokes equations, etc. In 2023, a benchmark[Hao et al., 2023], called PINNacle, is designed to assess the performance of 10 PINN variants, encompassed tasks derived from over 20 different PDEs spanning various domains, such as heat conduction, fluid dynamics, biology, and electromagnetic.

Previous research primarily benchmarks some popular methods such as DeepONet, FNO, U-Net, and PINN. However, a comprehensive evaluation of some other important NN methods is lacking. Additionally, different types of NN methods are often compared in the same setting without being classified and evaluated based on their distinct mathematical foundations. On the other hand, despite the selection of various types of PDE problems in previous work, some crucial equations, such as phase-field equations, have not been included in the evaluation tasks. Furthermore, challenging problems like PDEs with singular solutions have also been omitted. In summary, the development of the evaluation and comparison works for NN methods have not kept pace with the advancement of NN methods themselves.

In this work, we propose a comprehensive and systematic evaluation of various NN methods for solving diversified PDE problems, considering accuracy, efficiency and robustness.

Our evaluations are integrated into a benchmarks suite called PDENNEval, with some key features listed below.

• We provide a code reference that integrates 12 advanced NN methods, including six function learning-based NN methods: the Deep Ritz Method(DRM) [E and Yu, 2018], PINN [Raissi et al., 2019], Weak Adversarial Network (WAN) [Zang et al., 2020], Derivative-Free Loss Method (DFLM)[Han et al., 2020], Random Feature Method (RFM) [Chen et al., 2022], Deep Finite Volume Method (DFVM) [Cen and Zou, 2023], and six operator learning-based NN methods: U-Net[Ronneberger et al., 2015], Message Passing Neural PDE Solvers (MPNN)[Brandstetter et al., 2022], FNO [Li et al., 2020], DeepONet[Lu et al., 2021a], PhysicsInformed Neural Operator (PINO)[Li et al., 2021], Ushaped Neural Operators (U-NO)[Rahman et al., 2022]. We systematically evaluate the performance of these two categories in different setting.
• We carefully select 19 representative PDE problems from various fields such as fluid dynamics, materials science, finance, and electromagnetism to evaluate the performance of the aforementioned NN methods. These equations encompass some core mathematical challenges, including nonlinearity, singularity, high dimensionality and complex geometry. It's may worth mentioning that it is the first time NNmethods have been evaluated by using a PDE with a singular solution. More detailed descriptions for all PDEs are provided in Appendix A.
• We design a diverse and high-fidelity dataset, generated using high-precision conventional scientificcomputation methods. For instance, to generate the data for the phase-field equations, spatial discretization is accomplished using the spectral method, and time discretization is performed using the fourth-order RungeKutta method.
• We select a set of metrics that are better suited for each type of NN methods, including L2 norm relative error and maximum error, to assess accuracy, Lipschitz constant to evaluate robustness, and efficiency metrics encompassing training time, inference time, and convergence time.

The rest of the paper is organized as below. In Section 2, we outline the 12 representative NN methods and present 19 distinct PDE taskswhichwill be evaluated in PDENNEval. In Section 3, we present an overview of the datasets and metrics used in PDENNEval. In Section 4, we illustrate a standardized testing procedure, and present the results of the pertinent tests. Some brief concluding remarks are given in Section 5.

2 An Overview of Selected PDEs and NN Methods

In this section, we briefly introduce the PDE problems and NN methods integrated into the PDENNEval.

2.1 Selected PDE Tasks

A well-posed PDE system often involves three fundamental components: the equation itself, boundary conditions, and initial values. Let Ω ⊂ Rn be the spatial domain with boundary ∂Ω, and let [0, T] be the time interval. A general PDE system can be represented as:

where F represents a differential operator governing the behavior of the unknown function u, λ is an optional vector of parameters, appearing often as constant or variable coefficients of the differential operator, and f(x, t) is a vectorvalued function. This equation encapsulates the evolution of the system over the spatiotemporal domain. Boundary conditions are often expressed as:

where b is a given vector-valued function defined on ∂Ω. The boundary conditions B define how the solution behaves or is constrained at the boundary of the domain. At t = 0, the system's state is often specified by initial values:

where i is a given vector-valued function defined only on Ω.

Usually, we call it as a forward problem when seeking a function u to satisfy (1), (2) and (3), given the information Ω, λ, f, b, i. Conversely, we call it as an inverse problem if we aim to seek parameters λ, or the shape of Ω, or initial conditions, or other essential information about the system based on some observed data or outcomes. In this bench, since an inverse problem might be solved by training a NN using the same loss function as its corresponding forward problem, we will focus only on solving forward problems and encourage subsequent researchers to benchmark their respective inverse problems.

In our study, we aim to benchmark 19 forward problems of PDEs, representing pivotal mathematical models in physics and engineering that span various fields, including fluid dynamics, materials science, finance, and electromagnetism. Within the realm of fluid dynamics, our benchmark incorporates 11 equations: 1D Advection, 1D and 2D Burgers', 1D Diffusion-Reaction, 1D Diffusion-Sorption, 2D Darcy Flow, 2D Shallow-Water, 1D, 2D, and 3D Navier-Stokes, and 3D Euler equations. These equations model gas and liquid flow phenomena, finding applications across aerospace, weather prediction, oceanography, chemistry, and hydraulic engineering fields. Notably, while nine fluid-related PDEs have been evaluated in PDEBench, they were restricted to only three methods, whereas our assessment extends to twelve methods in total. In materials science, our selection includes the 1D and 2D Allen-Cahn equations and the 1D Cahn-Hilliard equation, employed in understanding material phase transitions, crystal growth, and surface evolution phenomena. Regarding finance, we encompass the 2D Black-Scholes-Barenblatt equation, crucial in the pricing of European options. In the domain of electromagnetism, our benchmark involves the Poisson equation and 3D Maxwell's equations, which model the behaviors of electromagnetic field and serve as fundamental elements in the realms of communication technology, electromagnetic device design, and optical systems. Specifically for the Poisson equation, we construct three distinct cases.

We select these PDEs to ensure that the benchmark introduces a wide range of mathematical challenges, including

• Nonlinear behavior: Numerous PDEs showcase nonlinear or even chaotic dynamics, wherein minor alterations in initial conditions result in significant divergence within the outcomes. Addressing these nonlinearities illuminates a method's capacity to capture intricate and sensitive dependencies within the system's behavior. In particular, our set of equations includes primarily nonlinear equations, such as the Navier-Stokes equations, phase field equations, and the Black-ScholesBarenblatt equation.
• High dimensionality: Various fields, including materials, finance, and quantum mechanics, present highdimensional PDEs. To evaluate the capability of NN methods in handling high-dimensional PDEs, we utilize high-dimensional (d > 3) Poisson equations to compare the performance of function learning-based NN methods.
• Singularity: Singularities in partial differential equations represent unconventional solution properties, often marked by unboundedness, discontinuities, or lack of differentiability at specific points. In our PDENNEval, we intentionally create a Poisson equation with a singular solution for assessment purposes.
• Complex geometry: PDEs of practical significance often extend beyond regular domains, underscoring the importance of evaluating NN methods' capability to handle irregular domain scenarios. In our study, focusing on the Poisson equation, we established a specific case involving the resolution of the equation over an Lshaped domain.

2.2 Selected NNMethods

An NN method heavily relies on its definition of the loss function L(θ) which depends on θ, the trainable parameters of the NN. Generally, a loss L(θ) for PDEs consolidates four essential components into a unified expression as follows:

where ωPDE, ωBC, ωIC, ωData are the weights of different loss terms, and a breakdown of each loss term is explained as below:

• PDE constraint loss (LPDE) enforces the PDE constraints within the domain Ω, defined as

where ̂ F is a variant of the original equations (1), including its strong, weak, stochastic forms; SPDE is the set of training points in the interior of Ω × [0, T], and N1 = #SPDE.

• Boundary condition loss (LBC) ensures adherence to specified boundary conditions at ∂Ω, defined as

where SBC is the set of training points sampled from the boundary ∂Ω × [0, T ], and N2 = #SBC.

• Initial condition loss (LIC) enforces accuracy at the initial time t = 0, defined as

where SIC is the set of training points sampled from the Ω, and N3 = #SIC.

• Data loss (LData) measures the discrepancy between predicted and observed data at designated points, defined as

where SData is the set of training points at which the value of the true solution u is known and N4 = #SData.

In general, physics-driven methods typically encompass only the first three components of the loss function, i.e., ωData = 0, while data-driven methods exclude the PDE constraint loss. However, in practice, certain physics-driven approaches encompass all loss components, blurring the distinction between data-driven and physics-driven paradigms. Hence, categorizing methods strictly into physics-driven and data-driven techniques becomes inadequate, given the strong coupling between these methodologies. The demarcation between these two approaches tends to be ambiguous.

In PDENNEval, we attempt to categorize NN methods for solving PDEs into two principal categories distinguished by their mathematical foundations. The first category employs elaborated NN to fit the solution function of the equation, rooted in the universal approximation theory[Hornik, 1991; Barron, 1994; He et al., 2022]. The second category utilizes NN to fit solution operators for a spectrum of PDEs, based on the universal approximation of linear or nonlinear operators by NN[Lu et al., 2021a]. Some details of these two methods are described as below.

Function Learning Based NN Methods

NNmethodologies for function learning largely fall under the domain of physics-driven methods, resembling traditional approaches in their goal of approximating solutions to PDEs. While traditional numerical methods often utilize piece-wise polynomials, NN methods leverage networks for this purpose. Typically, these methods involve training the networks using the four aforementioned loss components. However, similar to a traditional numerical method which solves a PDE

only using the equation itself, the boundary and initial conditions, a function learning type NN method often does not need to use exceptional data to train the NN. That is, ωData is often set to 0. Certainly, data might contribute to the boundary and initial value losses, but they are typically treated separately from the data loss component.

Function learning based methods differ mainly in their treatment of the PDE constraint loss item. Strong form-based methods like PINN and RFM integrate the PDE constraint loss utilizing the residual ̂ F = F - f. PINN and RFM both employ strong form-based loss functions, yet their optimization processes differ. PINN utilizes deep learning optimizers, while RFM, employing linear algebraic equations, fixes the parameters of the preceding layers to compute the final layer's parameters. Weak form-based methods encompass DRM, WAN, DFLM, and DFVM. DRM is rooted in the variational format, emphasizing the principle of minimal potential energy. WAN adopts the virtual-work-principle form, focusing on finding weak solutions to the equations. DFLM tailors its loss to match the equivalent stochastic-differential form of the original PDE. DFVM aligns its loss with the local-conservation form of the PDE.

Operator Learning Based NN Methods

Operator learning-based methods tackle solving a set of PDEs by learning some operators, categorized into autoregressive and non-autoregressive approaches. Autoregressive methods focus on learning mappings from solutions at preceding time steps to subsequent time steps, including FNO, MPNN, UNet, and U-NO. Non-autoregressive methods focus on learning mappings from auxiliary functions to solution functions, including DeepONet and PINO.

Autoregressive methods learn mappings from preceding time-step solutions to subsequent time-step solutions. These methods fundamentally simulate an image-to-image mapping, forecasting the future state based on the current state, which can be described as below :

These methods involve formulating the loss function to minimize the difference between the predicted solution at the next time step and the actual solution. This enables the network to learn the temporal evolution of the system. For instance, FNO focuses on parameterizing the integral kernel directly in Fourier space, enabling an expressive and efficient architecture for learning solution operators. MPNN operates on graph data by iteratively generating messages, aggregating neighbor information, and updating node representations. U-Net adapts an autoregressive approach for solving time-dependent PDEs. Its architecture predicts subsequent time step solutions based on preceding sequences of solutions, making it adept at capturing both local and global patterns for progressive solutions. U-NO employs a U-shaped architecture, mitigating memory constraints while efficiently mapping functions across different domains. By handling encoding and decoding with skip connections, it offers an effective solution for complex neural operator modeling.

Non-autoregressive methods involve learning an operator which maps the known information (λ, f, b, i) to solution functions u ∈ U, described as

In other words, these methods aim to learn mappings from right-hand side functions, boundary conditions, or other auxiliary information to the solution function. The loss function formulation for such mappings varies based on specific model designs and target functions, and may or may not involve the explicit PDE constraint loss term. For instance, DeepONet utilizes separate feed-forward neural networks (FNNs) to encode input functions and query locations, effectively mapping auxiliary functions directly to solutions. PINO blends data and physics information to learn resolution-invariant operators in Fourier space. It efficiently handles temporal PDEs by combining data-driven and physics-driven methodologies for accurate solutions across various resolutions and dimensions.

3 An Overview of Datasets andMetrics

In this section, we present two fundamental components of PDENNEval: the datasets and the metrics.

3.1 Datasets

The PDENNEval datasets consist of 16 data files, each representing a distinct PDE, excluding three instances of the Poisson equation which are solved using only the first two loss components in (4) (i.e. ωIC = ωData = 0). These 16 data files are named using the following abbreviations of PDEs.

In each data file, we store multiple samples, each of which is a solution of the related PDE. All solution samples are obtained by using traditional numerical methods to solve the PDE with different initial functions, sampled in a certain function space according to a certain distribution. Each sample contains the spatiotemporal coordinates of the underlying mesh and the values of the solution at those locations. Listed in Table 2 is a summary of all 16 data files. In this table, Nd is the number of spatial dimensions of the equation, and NS is the size of the data file. Moreover, we use Nsp, Nt, Nsa to denote the spatial resolution, the temporal resolution of the underlying mesh, and the number of samples, respectively. Note that, since the Darcy flow equation is time-independent, we have not provided the number Nt for the dataset DF2.

Some details about the traditional numerical methods for generating datasets are reported below. First, the following four datasets: AC1, AC2, CH1, and Bu2 are generated using the Matlab R2022a Software. More precisely, the samples in AC1, AC2 and CH1 are obtained by discretizing the spatial variable with a spectral method and discretizing the temporal variable with a fourth-order Runge-Kutta scheme. The samples in the dataset Bu2 are obtained by discretizing the spatial first-order derivatives with a first-order upwind scheme and discretizing the spatial second-order derivatives with a central difference scheme. Secondly, the datasets Ma3, BS2, Eu3 are generated using Python. To be more specific, the samples in Ma3 are obtained by discretizing the 3D Maxwell equation using the finite-difference time-domain method (FDTD, [Kunz and Luebbers, 1993]), and the resulting discrete algebraic system is solved with the leapfrog method. The BS2 dataset is sampled from some exact solutions of the BlackScholes-Barenblatt equation given in [Raissi, 2024]. The samples in Eu3 are obtained by solving the Euler equation utilizing the second-order HLLC scheme [Toro et al., 1994] for the inviscid part, the MUSCL schemes [Van Leer, 1979] for its inviscid counterpart, and a central difference scheme for the viscous part. The remaining nine datasets are related to PDEs in fluid dynamics, their samples are generated using the numerical methods which been used to generate the dataset in the PDEBench [Takamoto et al., 2022a].

All datasets are stored in HDF5 format, which is convenient to use and share and widely supported in a host of programs including Python, R and Matlab. The data and usage example can be found in https://github.com/Sysuzqs/PDENN Eval. More details of datasets can be found in Appendix C.

3.2 Metrics

In this section, we briefly introduce the metrics used in PDENNEval. To evaluate the accuracy, efficiency, and robustness of our selected NN methods, we collect 7 metrics which can be categorized into two groups: performance metrics and functional metrics. They are summarized in Table 3. Furthermore, we provide the definition of one metric from each of the aforementioned groups. The L2 norm relative error is defined as L2RE = ∥uθ - u∥2/∥u∥2. The Lipschitz constant Lips is defined as the smallest L satisfying

Here ∥·∥2 is the L2 norm, uθ is the predictive solution, and u is the reference solution. The definition of other metrics and their detailed description are given in Appendix D.

4 A Selection of Evaluation Results

In this section, we separately present partial experimental results for function learning-based NN methods and operator learning-based NN methods. Additional results can be found in Appendix F. For a fair comparison, all our experiments are conducted using the following system configuration: CPU : 2× Intel Xeon Gold 6230R @ 2.10GHz; GPU : 1 × NVIDIA TITAN RTX; Software : PyTorch@1.13.1, CUDA@11.6.

Recall that we have introduced 12 NN methods and 19 PDE tasks in Section 2. In our experiments, to test the performance of function learning NN methods, we introduce three additional PDE tasks PS, PL, PH which will be explained in details in Section 4.1. Apparently, in PDENNEval, not every method has been applied to all 19 PDE tasks. Listed in Table 4 is a summary of tasks coverage by our selected 12 NN methods.

From this table, we observe that PINN, DFVM, and RFM have been used to solve all 19 PDE problems; Meanwhile, the function learning methods WAN, DRM, and DFLM have only been applied to a part of these PDE problems. Following the vanilla methods, we use DRM [E and Yu, 2018] and DFLM [Han et al., 2020] only to solve symmetric and timeindependent problems PS, PL, PH and DF2. In the same way, we use WAN [Zang et al., 2020] only to solve scalar problems. In contrast, all operator learning type methods have not been used to solve the three Poisson tasks. Moreover, DeepONet, PINO and MPNN do not solve 3D problems, while the MPNN does not solve the time-independent problem DF2 either. It is worth mentioning that we do not mean that a specific NN method can not be applied to solve some specific PDE problems if it has not been done in our PDENNEval, it simply suggests that implementing this NN method for individual equation problems may pose challenges.

4.1 Results on Function LearningMethods

As already pointed out in Table 4, function learning-based NN methods can solve almost all kinds of PDE tasks. Here, we primarily present the experimental results of three tasks PH, PS and PL, which tackle three special mathematical challenges high-dimensionality, singularity, and complex geometry. The governing PDE of these three tasks is the Poisson equation

where Ω, f, g will be specified later.

1) High dimensional case (PH). In the first case, we evaluate the performance of function learning-based NN methods on high dimensional Poisson equations. Let Ω = [-1, 1]d , here d = 3, 5, 10, 20, 40, 80, 120 is the spatial dimension. The functions f and g are chosen so that u(x) = ( 1d ∑di=1 xi)2 + sin ( 1d ∑di=1 xi) , ∀x ∈ Rd , which is sufficiently smooth.

Listed in Table 5 are L2 norm relative error and training times of function learning methods for PH. From Table 5, we observe that all methods are immune to the curse of dimensionality. That is, the computation time of no method increases exponentially in terms of the dimension. Moreover, PINN exhibits superior performance in terms of the L2RE metric. The L2RE of the DFVM is almost as good as the PINN. In addition, on evaluation on the training time metric ttrain, the WAN excels for the 3, 5, 10, 20 dimensional cases, while DRM excels for the 40, 80 and 120 dimension cases. It is worth mentioning that since the RFM solution is obtained by solving a linear system of algebraic equations instead of training, no training time has been reported for this method.

2) Singular case (PS). In the second case, we consider the equation 11 with a singular solution. Precisely, let Ω = [0, 1]2 and f, g be selected to ensure u(x1, x2) = x21 if 0 ≤ x1 ≤ 12 and u(x1, x2) = (1 - x1)2 otherwise. Note that the exact solution u is continuous but not smooth at the location x1 = 1 2 .

Listed in Table 6 are L2RE and ttrain of function learning methods for this Poisson singular problem. We observe that in this case, the error of the DFLM is the smallest, followed by the DFVM. Moreover, the training time of the DFVM is the shortest. The error of the PINN and the RFM seems to be too big, which might be caused by the fact that they use the strong form PDE to construct their loss.

3) L-shape domain case (PL). In the third case, we choose Ω = [-1, 1]2 \(0, 1]2 and f = 1, g = 0. There is no analytical solution for this problem. We utilize the software COMSOL to generate a numerical solution whose configuration is shown in Fig1.

Listed in Table 7 are the L2 norm relative errors and training times for the problem PL. We observe that the error of the PINN is the smallest, followed by the DFVM. The training time of the DRM is the shortest.

Note that here, only results concerning L2RE and ttrain are reported, those concerning other metrics such as mERR, tinfer and tconv will be presented in the Appendix F.

4.2 Results on Operator LearningMethods

The six operator learning methods to be measured in PDENNEval are DeepONet, PINO, U-NO, FNO, U-Net, and MPNN, which are data-driven. Therefore, they are applied to solve 16 PDE problems except the three Poisson problems in the previous section, as shown in Table 4. Due to the space limitation, here we only present the evaluation results of the following three equations: the 1D Advection equation (Ad1), the 2D Allen-Cahn equation (AC2), and the 2D Darcy Flow equation (DF2). In our experiments, 90% of the data in each dataset is used for training, and 10% is used for testing. Moreover, here we only report the evaluation results in terms of the accuracy metric L2RE and three time metrics: the training time ttrain, the inference time tinfer, and the convergence time tconv.

Listed in Table 8 are the L2 norm relative error of the six operator learning methods. We observe that the U-NO achieves the best computational accuracy for problems Ad1 and AC2, while PINO achieves the best accuracy for the problem DF2. Note that we did not employ MPNN to solve the DF2.

Listed in Table 9 are the training, inference and convergence time. We observe that for all three PDE problems, DeepONet demonstrates noteworthy efficiency.

5 Concluding Remarks

In this study, we evaluate 6 function learning and 6 operation learning NN methods for solving 19 PDEs. Our experiments with function learning methods demonstrate that, for the PDEs with smooth solution, the PINN exhibits highest accuracy, and the DRM and WAN exhibits excellent training speed. In the case of PDE with non-smooth solution, the DFVM and WAN perform very well. Our experiments with operator learning methods reveal that, for autoregressive tasks, both U-NO and FNO perform admirably, while for non-autoregressive tasks, DeepONet exhibits the best performance.

Certainly, many other important NN methods, including variants of PINN, and many other fundamental PDE tasks such as multiscale problems in various fields have not been covered yet in PDENNEval. Our evaluation system is scalable and easily accessible, everyone is welcome to contribute data and/or codes to the https://github.com/Sysuzqs/PDENN Eval.

Developing PDE solvers with high accuracy, efficiency, and strong generalization is the pursuit of many computational scientists. The combination of traditional numerical methods and NN methods provides technical possibilities for achieving this pursuit. It can be anticipated that in the near future, there will likely be the emergence of many large NN models capable of efficiently solving lots of complex PDEs simultaneously (see e.g. [Yang et al., 2023a; Yang et al., 2023b] ). Therefore, in the future, we will pay much attention to the evaluation of large NN models for solving PDE problems.

Acknowledgments

This work was supported by the National Science and Technology Major Project (2022ZD0117805), by the National Natural Science Foundation of China under grants 92370113 and 12071496, and by the Natural Science Foundation of the Guangdong Province grant 2023A1515012097.

We thank graduate students Haolong Fan, Hongji Li and Changye He, all from Sun Yat-sen University, for their helps with numerical experiments.