Method for solving partial differential equations by physical information embedded graph neural network

Abstract

The application discloses a method for solving partial differential equations by embedding physical information into a graph neural network, comprising the following steps: S1, constructing a graph neural network GNN, replacing the physical relationship between the calculation nodes in the classical iterative solver with the edge message function of the GNN, and hard-coding the GNN components; S2, extracting local numerical state features; S3, constructing a compact parameter predictor omega-MLP; S4, embedding the solver iteration loop; S5, training driven by a pure physical loss; and constructing a pure physical loss function PDE residual. The application has the following advantages: physical information is embedded into the neural network structure in a hard constraint manner; the graph neural network structure directly corresponds to the difference format of the traditional iterative algorithm; the training process and the solving process are unified; the linear mixing mechanism provides a safe fallback capability; the compact MLP parameter quantity is irrelevant to the calculation grid scale; and the method is suitable for various discrete methods and physical problems.

Description

Technical Field

[0001] This invention relates to the field of computational methods, specifically to a method for solving partial differential equations by embedding physical information into a graph neural network. Background Technology

[0002] Partial differential equations (PDEs) are the core mathematical tools for describing the motion of continuous media in the physical world, forming the cornerstone of industrial simulation. Their solution techniques directly determine the design level and safety standards in key fields ranging from aerospace and automotive manufacturing to chip design and energy equipment. Specifically, computational fluid dynamics relies on the Navier-Stokes equations to optimize the aerodynamic shape of aircraft wings, high-speed train bodies, wind turbine blades, etc., and traditional simulations of a single operating condition often take tens of hours.

[0003] Each of these fields has spawned mature numerical simulation software, with an annual output value of tens of billions of dollars. Any breakthrough in the versatility and computational efficiency of partial differential equation solving technology will directly bring huge engineering value.

[0004] Currently, methods for solving partial differential equations for physical field simulation are mainly divided into two categories: traditional numerical methods and artificial intelligence-based methods. Traditional methods rely heavily on a set of manually set numerical control parameters. Furthermore, the same set of parameters may diverge or converge very slowly under different operating conditions, requiring engineers to undergo extensive trial and error to find suitable values. This manual parameter tuning method is not only inefficient but also severely limits the solver's ability to adaptively deploy across industrial scenarios. Therefore, various artificial intelligence methods have been developed to compensate for the shortcomings of traditional solvers. Among artificial intelligence-based methods, the closest existing technologies to this application are Physical Information Neural Networks (PINN), neural operator methods (FNO / DeepONet), and purely data-driven graph neural network surrogate models. All of these methods have their limitations.

[0005] Physical Information Neural Networks (PINs) are a method that uses fully connected neural networks to solve Partial Differential Equations (PDEs) by using the residuals of PDEs as a loss function. The workflow is as follows: the neural network receives spatial coordinates (x, y) as input and outputs predicted physical quantities at that point (e.g., pressure P, velocity u); then, it automatically differentiates and calculates the partial derivatives of the network output with respect to the input to construct the PDE residuals; finally, these residuals, along with the boundary condition errors, are used as the loss function to train the network. However, this method has the following drawbacks:

[0006] (1) The physical equation is a soft constraint. The physical equation is added to the loss function in the form of a penalty term. The network may not be able to accurately satisfy the equation during training, especially when dealing with strongly nonlinear problems or regions with large gradients, the training difficulty increases significantly.

[0007] (2) Retraining is required when the operating conditions change. When the operating conditions change, the network usually needs to be retrained from scratch.

[0008] (3) The reasoning method is point-by-point calculation. It is necessary to sample coordinate points in the computational domain and reason one by one to obtain the overall solution.

[0009] (4) No convergence guarantee. Training may not converge or may converge to an incorrect solution, which the user cannot know in advance.

[0010] The neural operator method uses a neural network to learn the mapping from the parameter space to the solution space (NN = operator), that is, to learn a function G: input parameters → PDE solution. Its workflow is as follows: first, a large number of simulations under different conditions are run using a traditional numerical solver to generate a massive amount of "parameter-solution" pairing data as a training set; then, a neural network model is trained to minimize the error between its predicted solution and the simulation label. After training, the inference speed is very fast. However, it has the following shortcomings:

[0011] (1) A large amount of training data is required. Usually, a high-fidelity reference solution needs to be solved in advance as training data. For complex problems, the computational cost of generating training data is high.

[0012] (2) The number of model parameters is usually large.

[0013] (3) The prediction results may violate basic physical laws such as conservation of mass and conservation of momentum.

[0014] Data-driven graph network proxy models, represented by MeshGraphNets, are based on the core idea of ​​mapping a physical system's mesh to a graph structure and then using a graph neural network (GNN) to learn the mapping from the current state to the next time step. The workflow is as follows: First, a large number of simulations under different operating conditions are run using traditional numerical solvers to generate massive amounts of high-fidelity simulation data (i.e., state sequences) as training labels; then, a GNN model is trained to minimize the error between its predicted next time step state and the simulation labels. While this method offers fast inference speed, it suffers from the following drawbacks:

[0015] (1) Usually, a large number of high-precision traditional simulations are required to generate training data, which consumes a lot of computing resources and time.

[0016] (2) Limited generalization ability. When encountering new working conditions with a different distribution than the training data, the model prediction may deviate from the true physical laws.

[0017] (3) When the model makes a wrong prediction, there is no mechanism to fall back to the classical solver.

[0018] In summary, all the above schemes attempt to "replace" a certain function of the classical solver (representing solutions, operators, iteration steps, or simulators) with neural networks. This "replacement" strategy may lead to: The convergence guarantees and error estimation theories accumulated by classical methods were not preserved; Neural networks lack a safe rollback mechanism when errors occur; The number of model parameters is usually large and increases with the problem size; The versatility across PDE types and operating conditions needs to be improved.

[0019] To address the aforementioned problems, this invention presents a method for solving partial differential equations by embedding physical information into a graph neural network. Summary of the Invention

[0020] To address the shortcomings of existing technologies, this invention provides a method for solving partial differential equations by embedding physical information into graph neural networks. Its core lies in hard-coding the discrete format of physical governing equations, such as the Reynolds equation and the Navier-Stokes equation, into the message passing process of the graph neural network. Combined with a purely physical-driven training and online fine-tuning mechanism, this fundamentally overcomes the deficiencies of other types of neural networks in solving partial differential equations. Furthermore, it strictly maps the computational structure of traditional numerical iterative solvers onto the computational graph of the graph neural network to solve the aforementioned problems.

[0021] To achieve the above objectives, this invention provides the following technical solution: a method for solving partial differential equations using a graph neural network embedded with physical information. This neural network not only retains the computational structure of traditional numerical iterative solvers but also learns the numerical control parameters of classical solvers, whereas the numerical control parameters of traditional algorithms typically rely on manual adjustment. The iterative format of a traditional iterative solver is as follows: ,in These are numerical control parameters. This invention will fix the parameters. Replace with a learning function. ,in This is a local numerical state feature vector. This method utilizes a compact multilayer perceptron (MLP) with a parameter count independent of the computational grid size. This paradigm shift, from learning the mapping from input to output to learning the solution process of traditional algorithms, overcomes the shortcomings of existing technologies. The specific steps of this method are as follows: Step S1: Construct a graph neural network (GNN), and directly replace the physical relationships (flux terms in the discrete format) between computation nodes in the classic iterative solver with the edge message functions of the GNN, which are hard-coded into the components of the GNN. The mapping relationship between the traditional algorithm and the GNN is as follows: (1) Map all computation nodes in the discrete format to nodes of the GNN, and each node carries physical field values, residuals, local physical coefficients and spatial coordinates as node features; (2) Establish directed edges of GNN between all pairs of nodes that have direct interaction in the discrete format, and use the surface average physical coefficients carried by the edges as edge features. (3) The inter-node flux formulas given in the discrete format (such as Fourier heat flux in heat conduction, momentum flux in fluid mechanics, elastic force flux in structural mechanics, etc.) are directly hard-coded into the side message function m(j→i) of GNN — that is, the function form of the side message is determined by the discrete format of the partial differential equation and does not contain any trainable parameters. (4) The flux divergence calculation in the discrete format (corresponding to the conservation laws of mass, momentum or energy) is hard-coded into the node aggregation function of GNN, that is, the sum of incoming edge messages minus the sum of outgoing edge messages; (5) The iterative update formula in the discrete scheme is retained as the node update function of the GNN, but the original fixed relaxation factor omega is replaced with the adaptive value omega_MLP(i) predicted by the MLP; this is the part of the entire GNN containing trainable parameters. Through the above construction, the resulting GNN is itself a complete partial differential equation iterative solver. Each layer of message passing in the GNN is equivalent to one iteration of the classical solver under the discrete scheme. The graph structure, edge message function, and node aggregation function of the GNN are all explicitly determined by the selected discrete scheme and do not need to be learned through data training.

[0022] Step S2: Extract local numerical state features. During the solver iteration process, a local numerical state feature vector s(x) is extracted for each computational grid point. This feature vector includes: the physical field values ​​at that point (e.g., pressure P, film thickness H), the partial differential equation residual R, the absolute value of the residual |R|, the sum of local diffusion coefficients, and the spatial coordinates (X, Y). The feature design principle is to include all physical quantities affecting local convergence stability, enabling the MLP to determine the convergence risk and optimal step size for that point based solely on the local numerical state.

[0023] Step S3: Construct a compact parameter predictor (omega-MLP). Construct a compact multilayer perceptron as the parameter predictor N_phi: s(x)→theta(x). Its preferred architecture is: Input layer: 7-8 neurons (corresponding to the dimension of the feature vector); Hidden layer: 12-16 neurons, activation function is ReLU; Output layer: 1 neuron, activation function is tanh saturation function; The total number of parameters |phi| is between 100 and 300, preferably 8→16→1, for a total of 161 parameters; key property: the number of parameters in this MLP is independent of the computational grid size, and the same MLP can serve 17^2 grids, 129^2 grids, and even grids of arbitrary size. The tanh saturation function of the output layer ensures that the relaxation factor is restricted within the stability domain, with omega belonging to (0, omega_max), and omega_max=0.6 being a hard upper bound for stability.

[0024] Step S4: Embedded Solver Iteration Loop. In each iteration of the solver, the adaptive numerical control parameters for each grid point are predicted using the MLP, replacing the original fixed parameter values. For each grid point (i,j), the iterative update format is: The rest of the solver remains completely unchanged: the partial differential equation discretization scheme, boundary condition handling, convergence criteria, multi-grid hierarchy (FAS V-cycle / W-cycle), and FFT to accelerate elastic deformation calculations. This "embedded rather than replaced" design ensures that all convergence theories of classical solvers still apply; MLP only affects the convergence speed, not the correctness of the convergence results (correctness decoupling).

[0025] Step S5: Pure physical loss-driven training (zero data dependency); construct a pure physical loss function (PDE residual) without relying on any pre-computed reference solution data.

[0026] Global load balance constraint gradients can be calculated using either the finite difference method (SPSA) or automatic differentiation, with the Adam optimizer used to update MLP parameters. Training requires only PDE residuals and no labeled data.

[0027] Step S6: Three optional deployment modes: Mode A (Iterative Acceleration): Predicts pointwise adaptive omega(x) in the standard FAS V-cycle, replacing the fixed relaxation factor, resulting in a speedup of 1.2-5 times, with complete physical guarantees.

[0028] Mode B (One-Step Inference): Freezes the trained MLP and directly outputs the converged solution through forward propagation of a fixed-layer GNN, without the need for iteration and convergence judgment. It achieves a speedup of approximately 58 times on a 129^2 grid (compared to approximately 29.2 seconds for traditional FAS).

[0029] Mode C (Online Fine-tuning with Few Samples): Fine-tunes only 100-300 parameters of the MLP (grid / map / kernel freeze), with 10-50 steps of gradient update (2-10 seconds), enabling rapid adaptation to new operating conditions.

[0030] Furthermore, graph neural network physical flux hard coding: Step A1: Graph Mapping. Construct a directed bipartite graph G=(V,E) on an N×N finite difference grid. The node set V corresponds to N^2 grid points. Each internal node has 4 incoming edges and 4 outgoing edges. The edge characteristic is the surface-average diffusion coefficient.

[0031] Step A2: Edge messages are hardcoded as physical flux. The message for edge (i→j) is defined as: This message is mathematically equivalent to the Poiseuille flux term (the flux component of the Reynolds equation). This mapping is an exact algebraic identity in a five-point finite difference discretization scheme, not an approximation—the equivalence holds independently of the MLP parameters at machine precision.

[0032] Step A3: Node aggregation is calculated as flux divergence. The message aggregation value of node v = sum(incoming edge messages) - sum(outgoing edge messages), which is mathematically equivalent to flux divergence calculation and corresponds exactly to the Poiseuille part in the Reynolds residual.

[0033] Furthermore, the Fallback safety rollback mechanism: Set the linear mixing ratio f to [0,1]: omega_eff(x) = omega_fixed + f * (omega_MLP(x) - omega_fixed) When f=0: The system is completely equivalent to a traditional fixed-parameter solver, with a numerical error of less than 10^-12 (machine accuracy). When f=1: It fully utilizes the adaptive parameters predicted by the MLP, achieving maximum speedup. When f is between 0 and 1: It smoothly transitions between a safe baseline and maximum speedup. This mechanism provides industrial-grade safety: the solver can always operate at a safe baseline of f=0, gradually increasing f as trust increases. This is a capability that existing NN methods such as PINN and FNO do not possess.

[0034] Furthermore, multi-parameter joint learning: Expanded to a combination of three MLPs, with a total of 279 parameters: omega-MLP (161 parameters, 8→16→1): Predicts node-wise adaptive relaxation factors; eta-MLP (65 parameters, 4→8→1): Predicts load regulation rate from global features; gamma-MLP (53 parameters, 4→8→1): Predict coarse-grid corrected damping from residual statistics; The three MLPs are jointly trained within the FAS multigrid framework, and the total loss is PDE residual + load constraint + damping regularization.

[0035] This invention offers the following advantages: physical information is embedded into the neural network structure in a hard-constraint manner; unlike PINN, which uses physical equations as soft constraints in the loss function, this invention directly hard-codes the flux formula, explicitly given in the discrete form of the partial differential equation, into the side message function of the graph neural network (GNN). The calculation of physical flux becomes part of the network structure rather than the training objective. This ensures that each message pass in the GNN executes an exact algebraic identity in the discrete form, rather than a data-driven approximation.

[0036] The graph neural network (GNN) structure directly corresponds to the difference scheme of traditional iterative algorithms. The GNN graph structure of this invention is not learned from data, but rather a direct mapping of the discrete scheme of traditional solvers: the difference relationship determines whether there are edges connecting nodes, and the flux calculation formula determines the message function on the edges. The message passing process of GNN has the same mathematical structure as the iterative calculation process of traditional iterative solvers.

[0037] The training and solution processes are unified. Unlike existing methods that employ a two-stage approach of training the model first and then inferring and solving, this invention alternates between parameter updates and iterative solutions of partial differential equations in the same loop for the compact MLP. Each iteration of the solver provides the MLP with training signals (PDE residuals), and the parameter updates of the MLP provide better relaxation factors for the next iteration. The MLP is also trained when the solution converges, eliminating the need for pre-prepared training datasets.

[0038] A safe fallback capability is provided through a linear hybridization mechanism. By setting the hybridization ratio f∈[0,1], the solver can continuously transition between a traditional fixed-parameter mode (f=0) and a fully adaptive mode (f=1). When f=0, the solver degenerates into a traditional iterative solver, and the correctness of the solution is not affected by the MLP parameters.

[0039] The number of parameters in a compact MLP is independent of the computational grid size. The same MLP model (e.g., 8→16→1, a total of 161 parameters) can serve grids of different sizes, and the number of parameters does not increase with grid refinement.

[0040] It is applicable to a variety of discretization methods and physical problems. The mapping relationship between GNN and partial differential equation difference schemes (computation node → graph node, flux → edge message, divergence → node aggregation, iterative update → node update) is applicable to various discretization methods such as finite difference method, finite volume method, finite element method and meshless particle method, as well as physical problems with flux-form discretization schemes such as heat conduction, fluid mechanics, and structural mechanics. Attached Figure Description

[0041] Figure 1 This invention uses elastohydrodynamic lubrication as an example to calculate the neural network training and calculation process flowchart.

[0042] Figure 2 The flowchart shows the traditional iterative algorithm for loop calculation, taking elastohydrodynamic lubrication as an example. Detailed Implementation

[0043] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0044] Example 1: Solving point contact problems in elastohydrodynamic lubrication (EHL) (EHL, elliptic nonlinear PDE) Please see Figure 1 This embodiment demonstrates the application of the present invention to coupled multiphysics nonlinear problems. The governing equations for the steel ball-plate point contact EHL problem are the coupling of the Reynolds equation (fluid lubrication), the Boussinesq elastic integral (solid deformation), and the load equilibrium condition (force equilibrium). Physical parameters: steel ball radius R = 0.02 m, load W = 20 N, and entrainment velocity Us = 1.0 m / s.

[0045] Mapping Relationship: On a 129×129 finite difference grid, grid points are used as GNN nodes, and bidirectional edges are established between 4 adjacent grid points. Edge messages are hard-coded as Poiseuille flux (directly given by the discrete format of the Reynolds equation). Node aggregation is the sum of the four incoming edge messages minus the sum of the four outgoing edge messages, corresponding to the Poiseuille flux divergence. The source term is the Couette flux. Node updates use the Newton-Raphson format, where the relaxation factor is predicted by omega-MLP. omega-MLP takes 8-dimensional features [P, H, R, |R|, eps_sum_x, eps_sum_y, X, Y] as input and outputs a pointwise adaptive relaxation factor. Training is conducted online (training and solving are synchronized), and the loss function is the Reynolds residual mean square value + load balance constraint.

[0046] Results: Correct solutions with load errors below 1% were obtained under six different operating conditions (load W = 10-50 N, entrainment speed Us = 0.5-5.0 m / s). The learned omega distribution exhibits a physically interpretable spatial structure: it automatically decreases in the pressure peak region to ensure stability in the nonlinear region, and increases in the smooth inlet region to accelerate convergence. Fallback validation: When f = 0, the values ​​are completely consistent with those of traditional solvers (error < 10^-12).

[0047] Example 2: Steady-state heat conduction (elliptic linear PDE) The governing equation is ∇·(k(x,y)∇T)=q(x,y), where k(x,y) is the non-uniform thermal conductivity (e.g., containing Gaussian hotspots), and q(x,y) is the distributed heat source.

[0048] Mapping Relationship: On an N×N finite difference grid, grid points are used as GNN nodes, and bidirectional edges are established between 4 adjacent grid points. The edge message is hard-coded as the Fourier heat flux m(j→i)=k_ij·(T_j-T_i) / Δ^2 (where k_ij is the surface average thermal conductivity, which is directly given by the discrete format of the heat conduction equation). Nodes are aggregated into four-sided heat flux divergence, with the source term S_i=q_i. Node updates use the Jacobi or SOR format, and the relaxation factor is predicted by omega-MLP. The omega-MLP takes 7-dimensional features [T, R, |R|, k_sum_x, k_sum_y, X, Y] as input and outputs a point-by-point adaptive relaxation factor. Since the heat conduction equation is a linear PDE, the optimal relaxation factor is close to a constant. The b2 bias of the MLP is calibrated to 8.0 to make the initial output close to the theoretical optimum.

[0049] The mapping relationships in this embodiment correspond completely to those in Embodiment 1: physical quantity P → T, flux coefficient eps → k, source term Couette flux → heat source q. The structures of the four mapping relationships are exactly the same; only the physical meaning of the physical quantities and the calculation formulas of the coefficients change with the governing equations.

[0050] Example 3: Cavity-driven flow (Navier-Stokes equations, hybrid PDE) The governing equations are the incompressible Navier-Stokes equations, discretized using the stream function-vorticity (ψ-ω) scheme. The wall vorticity is represented by the Thom formula, and the convection term is represented by the upwind mixing scheme.

[0051] Mapping Relationship: On a uniform grid, grid points are used as GNN nodes, and bidirectional edges are established between adjacent grid points. Edge messages are hard-coded as convection-diffusion fluxes, including convection terms (given directly from the discrete scheme by the product of the stream function gradient and vorticity) and diffusion terms (vorticity Laplacian divided by the Reynolds number). Node aggregation is the convection-diffusion flux divergence, corresponding to the vorticity transport equation. Node updates adopt an explicit time-progression scheme, and its upwind mixing factor q∈[0,1] is learned by omega-MLP (here, MLP learns the upwind weights rather than the relaxation factor, reflecting the flexibility of PI-GNN in learning different numerical parameters).

[0052] Results: When Re=100, the MLP learns that q=0 (pure central difference) is optimal, with a central linear velocity error of only 1.14% (compared to the Ghia benchmark).

[0053] Example 4: Cylindrical Turbidity (Navier-Stokes Equations, Convection-Dominated PDE) The governing equations are the same as in Example 3, but the domain contains a cylindrical obstacle (diameter D=0.1), which is handled using the Brinkman volume penalty immersion boundary method.

[0054] Mapping Relationships: In the graph structure, the nodes and edges of the fluid region are the same as in Example 3. Nodes inside the cylinder are labeled as obstacle types. In addition to the regular convection-diffusion flux, the edge messages of nodes near the cylinder surface also include an immersion boundary penalty flux term (given directly from the discrete format by the penalty coefficient × local velocity difference). omega-MLP learns the optimal upwind mixing factor q, and the input features also include the distance from the cylinder surface. The remaining mapping relationships are the same as in Example 3.

[0055] Example 5: Multi-parameter Joint Learning. Building upon Example 1, the parameter predictor is expanded to a combination of three independent MLPs: omega-MLP (8→16→1, 161 parameters), eta-MLP (4→8→1, 65 parameters, learning load regulation rate from global features), and gamma-MLP (4→8→1, 53 parameters, learning coarse-grid corrected damping from residual statistics). The three MLPs have a total of 279 parameters. The mapping relationships, graph structure, edge messages, and node aggregation are completely identical to those in Example 1; only the learnable parameters in node updates are expanded from one (relaxation factor) to three (relaxation factor + load regulation rate + coarse-grid corrected damping).

[0056] Summary of Examples: The above examples cover four types of PDEs: elliptic nonlinear (EHL), elliptic linear (heat conduction), hybrid (cubic cavity flow), and convection-dominated (flow around a cylinder), involving three physical fields: lubrication mechanics, thermodynamics, and fluid mechanics. The core mapping structure of all examples is completely consistent—computation node → GNN node, flux formula → edge message (hard-coded), flux aggregation → node aggregation, iterative update → node update (parameters learned by MLP)—only the physical quantities and coefficient formulas change with the governing equations.

[0057] In addition to using a multilayer perceptron (MLP) as the parameter predictor, the following neural network architectures can also be used: (a) KAN (Kolmogorov-Arnold Networks): KAN uses B-spline basis functions as activation functions, and the learned B-spline coefficients have direct physical interpretability—the B-spline coefficients in different intervals correspond to the relaxation factor values ​​in different physical parameter intervals. KAN's spline functions naturally support extrapolation, and may have better stability when extrapolating to conditions not covered by training.

[0058] (b) Transformer: It uses a self-attention mechanism as a parameter predictor and automatically captures long-range dependencies between computing nodes through attention weights, thereby enhancing its adaptability to problems with nonlocal physical effects (such as elastic deformation and electromagnetic fields).

[0059] (c) Convolutional Neural Network (CNN): For problems on regular grids, small CNNs (such as 3×3 or 5×5 convolution kernels) can be used instead of MLPs to reduce the number of parameters by taking advantage of the translational equivariance of convolution.

[0060] (d) Gaussian Process: Using a Gaussian process as a parameter predictor, it provides an estimate of the uncertainty of the prediction while providing parameter prediction, which can be used to adaptively control the mixing ratio f.

[0061] 2. Alternative implementations of the discrete method.

[0062] The mapping relationship between traditional algorithms and neural networks described in this invention is not limited to the finite difference method. The corresponding mapping methods of the following discrete methods are all within the scope of protection of this invention:

[0063] (a) Finite Volume Method (FVM): The center of the control volume is used as the graph node, and graph edges are established between adjacent control volumes sharing a surface. The edge message function is determined by the numerical flux formula given by the selected Riemann solver (such as Roe format or HLLC format). Node aggregation corresponds to flux balance on the control volume.

[0064] (b) Finite Element Method (FEM): Graph nodes are set up with element nodes or integration points as graph nodes and graph edges are established between nodes of shared elements. The edge message function is determined by the corresponding elements of the element stiffness matrix. Node aggregation corresponds to the assembly of the global stiffness matrix.

[0065] (c) Meshless particle method (such as SPH, MPS, etc.): Particles are used as graph nodes, and graph edges are established between particle pairs in the support domain. The edge message function is determined by the selected kernel function (such as Wendland kernel, Cubic Spline kernel) and its derivative. Node aggregation corresponds to the summation of particle interactions in the SPH format.

[0066] (d) Spectral method: Using collocation points or modal coefficients as graph nodes, graph edges are established based on coupling relationships in the spectral space.

[0067] 3. Numerical control parameters other than the relaxation factor.

[0068] Although this invention primarily uses the relaxation factor omega as its embodiment, the PI-GNN paradigm can learn any numerical control parameter in a classical iterative solver, including but not limited to:

[0069] (a) Multi-grid cycle type parameters: number of V-cycle or W-cycle layers, number of pre-smoothing and post-smoothing cycles per layer.

[0070] (b) Coarse mesh correction damping coefficient: controls the confidence level of coarse mesh correction.

[0071] (c) Adjustment step size of load regulation rate or global constraint.

[0072] (d) Artificial viscosity coefficient: the intensity of artificial dissipation in shock wave capture or convection-dominated problems.

[0073] (e) Time step factor: Adaptive time step selection in transient solutions.

[0074] (f) Multi-parameter joint learning: Use multiple independent compact MLPs to learn different control parameters, or unify all parameters by outputting a single MLP (the output layer consists of multiple neurons).

[0075] 4. Alternative implementation of training methods.

[0076] (a) Online training (training and solving simultaneously): As described in the examples, the MLP parameters are updated in real time during the iterative solving process. This is suitable for scenarios where each solving task is optimized independently.

[0077] (b) Offline pre-training + online fine-tuning: First, pre-train the MLP on a representative set of work cases, and then use the pre-trained weights as the initial values ​​for new work cases. Only a small amount of online updates are needed for adaptation. This is suitable for scenarios with a large number of similar problem-solving tasks.

[0078] (c) Meta-Learning: Train the MLP on multiple operating conditions to enable it to learn how to quickly adapt to new operating conditions, so that it can reach the optimum with only a few gradient updates on new operating conditions.

[0079] (d) Evolution Strategies: These strategies replace SPSA for gradient estimation. They estimate the gradient by sampling multiple random perturbation directions, which can reduce gradient noise and are suitable for scenarios that require higher precision gradient estimation.

[0080] (e) Bayesian optimization: Instead of gradient-based optimization methods, Bayesian optimization searches for the optimal MLP weights in the parameter space, which is suitable for scenarios where gradients are difficult to obtain or where there is a lot of noise.

[0081] 5. Alternative construction of the loss function.

[0082] (a) Using only local PDE residuals: For linear PDEs that are automatically satisfied by global conservation (such as the Poisson equation and the heat conduction equation), only the mean square value of the local residuals can be used as the loss.

[0083] (b) Add boundary condition loss term: explicitly add the violation of boundary conditions to the loss function to accelerate convergence near the boundary.

[0084] (c) Add historical information loss term: Use the residual decay trend of the previous few iterations as an additional training signal.

[0085] (d) Using only global conservation loss: For problems dominated by global behavior, only the conservation error of global quantities (such as total force, total energy) can be used as the loss.

[0086] 6. Alternatives to embedding physical information into GNNs.

[0087] (a) Incomplete hard coding: Some parameters in the side message function (such as the empirical constant in the flux coefficient) can be set as trainable parameters and learned together with the relaxation factor.

[0088] (b) Multi-level cascaded graph: A graph is constructed at each level of the FAS multigrid. The edge message function of the coarse graph uses the same mapping rule (determined by the discrete format of that level). Different levels are connected by restrictive and extended edges.

[0089] (c) Heterogeneous edge types: Different edge message functions are used according to different types of physical flux (convection, diffusion, pressure, etc.), and they are summed separately and then combined when the nodes are aggregated.

[0090] 7. Alternative implementations of deployment methods.

[0091] (a) The trained MLP model can be exported as a standalone file (such as ONNX format) and embedded into other simulation software or hardware without relying on the original training environment.

[0092] (b) Integration of multiple MLP models: train multiple MLPs with different initializations and average or vote on their predictions during inference to improve robustness.

[0093] (c) Layered deployment: Use the relaxation factor predicted by MLP in critical areas (such as pressure peaks and boundary layers) and use the fixed relaxation factor in non-critical areas to balance computational overhead and acceleration effect.

[0094] The preferred embodiments of the present invention disclosed above are merely illustrative of the invention. These preferred embodiments do not exhaustively describe all details, nor do they limit the invention to the specific implementations described. Clearly, many modifications and variations can be made based on the content of this specification. This specification selects and specifically describes these embodiments to better explain the principles and practical applications of the invention, thereby enabling those skilled in the art to better understand and utilize the invention. The invention is limited only by the claims and their full scope and equivalents.

Claims

1. A method for solving partial differential equations by embedding physical information into a graph neural network, characterized in that: The iterative format of traditional iterative solvers Fixed parameters Replace with learning function ,in This is a local numerical state feature vector. It is a compact multilayer perceptron (MLP) with a parameter count independent of the computation grid size; Specifically, the following steps are included: Step S1: Construct a graph neural network (GNN), directly replacing the physical relationships between computation nodes in the classic iterative solver with the edge message functions of the GNN, and hard-coding them into the components of the GNN; Step S2: Extract local numerical state features; During the solver iteration process, extract the local numerical state feature vector s(x) for each computation grid point; Step S3: Construct a compact parameter predictor omega-MLP; construct a compact multilayer perceptron as the parameter predictor N_phi: s(x)→theta(x); Step S4: Embed the solver iterative loop; in each iteration of the solver, use MLP to predict the adaptive numerical control parameters for each grid point and replace the original fixed parameter values; for each grid point (i,j), the iterative update format is: The rest of the solver remains completely unchanged: partial differential equation discretization scheme, boundary condition handling, convergence criteria, multi-grid hierarchy (FAS V-cycle / W-cycle), and FFT to accelerate elastic deformation calculation; Step S5: Pure physical loss-driven training; construct the pure physical loss function PDE residuals, without relying on any pre-computed reference solution data; Step S6: Three optional deployment modes: Iterative acceleration: Predict pointwise adaptive omega(x) in the standard FAS V-cycle, replacing the fixed relaxation factor, resulting in a speedup of 1.2-5 times, with complete physical guarantees; One-step reasoning: Freeze the trained MLP and directly output the converged solution through forward propagation of a GNN with a fixed number of layers, without the need for iteration and convergence judgment; Online fine-tuning with few samples: Only 100-300 parameters of the MLP are fine-tuned, and gradient updates are performed in 10-50 steps to achieve rapid adaptation to new operating conditions.

2. The method for solving partial differential equations by embedding physical information into a graph neural network according to claim 1, characterized in that: The mapping relationship between traditional algorithms and GNNs is as follows: (1) Map all computation nodes in the discrete format to nodes of GNN, with each node carrying physical field values, residuals, local physical coefficients and spatial coordinates as node features; (2) Establish directed edges of GNN between all pairs of nodes that have direct interaction in the discrete format, and use the surface average physical coefficients carried by the edges as edge features. (3) The inter-node flux formula given in the discrete format is directly hard-coded into the side message function m(j→i) of GNN. That is, the functional form of the side message is determined by the discrete format of the partial differential equation and does not contain any trainable parameters. (4) The flux divergence calculation in the discrete format is hard-coded into the node aggregation function of GNN, that is, the sum of incoming edge messages minus the sum of outgoing edge messages; (5) The iterative update formula in the discrete scheme is retained as the node update function of GNN, but the original fixed relaxation factor omega is replaced with the adaptive value omega_MLP(i) predicted by MLP. Through the above construction, the obtained GNN is itself a complete partial differential equation iterative solver. Each layer of message passing in GNN is equivalent to one iteration of the classical solver under the discrete scheme. The graph structure, edge message function and node aggregation function of GNN are all explicitly determined by the selected discrete scheme and do not need to be learned through data training.

3. The method for solving partial differential equations by embedding physical information into a graph neural network according to claim 2, characterized in that: In step S2, the feature vector includes: the physical field value of the point, the partial differential equation residual R, the absolute value of the residual |R|, the sum of the local diffusion coefficients, and the spatial coordinates (X,Y). The feature design principle is to include all physical quantities that affect the local convergence stability, so that the MLP can determine the convergence risk and optimal step size of the point based solely on the local numerical state.

4. The method for solving partial differential equations by embedding physical information into a graph neural network according to claim 3, characterized in that: In step S3, the architecture of the parameter predictor is as follows: Input layer: 7-8 neurons; Hidden layer: 12-16 neurons, activation function is ReLU; Output layer: 1 neuron, activation function is tanh saturation function; The total number of parameters |phi| is between 100 and 300; the number of parameters of this MLP is independent of the computation grid size, and the same MLP can serve 17^2 grids, 129^2 grids, and even grids of arbitrary size; the tanh saturation function of the output layer ensures that the relaxation factor is restricted within the stability domain, omega belongs to (0, omega_max), and omega_max=0.6 is a hard upper bound for stability.

5. The method for solving partial differential equations by embedding physical information into a graph neural network according to claim 4, characterized in that: In step S5, the global load balance constraint gradient can be calculated using the finite difference method SPSA or automatic differentiation, and the Adam optimizer is used to update the MLP parameters; training only requires PDE residuals and no labeled data is needed.

6. The method for solving partial differential equations by embedding physical information into a graph neural network according to claim 5, characterized in that: Hard-coding of physical flux in graph neural networks: Step A1: Graph mapping; Construct a directed bipartite graph G=(V,E) on an N×N finite difference grid; The node set V corresponds to N^2 grid points; Each internal node has 4 incoming edges and 4 outgoing edges; The edge characteristic is the surface average diffusion coefficient; Step A2: Edge messages are hardcoded as physical flux; the message definition for edge (i→j) is: This message is mathematically equivalent to the Poiseuille flux term; this mapping is the exact algebraic identity of the five-point finite difference discrete scheme. Step A3: Node aggregation is the flux divergence; the message aggregation value of node v = sum(incoming edge messages) - sum(outgoing edge messages), which is mathematically equivalent to the flux divergence calculation, corresponding to the Poiseuille part in the Reynolds residual.

7. The method for solving partial differential equations by embedding physical information into a graph neural network according to claim 6, characterized in that: A fallback mechanism is used: the linear blending ratio f is set to [0,1]: omega_eff(x) = omega_fixed + f * (omega_MLP(x) - omega_fixed) When f=0: the system is completely equivalent to a traditional fixed-parameter solver with a numerical error of less than 10^-12; when f=1: it fully utilizes the adaptive parameters predicted by the MLP to achieve maximum speedup; when f is between 0 and 1: it provides a smooth transition between the safety baseline and maximum speedup; this mechanism provides industrial-grade safety assurance: the solver can always operate at the safety baseline of f=0, and gradually increase f as trust increases.

8. The method for solving partial differential equations by embedding physical information into a graph neural network according to claim 7, characterized in that: Multi-parameter joint learning is employed: Expanded to a combination of three MLPs, with a total of 279 parameters: omega-MLP (161 parameters, 8→16→1): Predicts node-wise adaptive relaxation factors; eta-MLP (65 parameters, 4→8→1): Predicts load regulation rate from global features; gamma-MLP (53 parameters, 4→8→1): Predict coarse-grid corrected damping from residual statistics; The three MLPs are jointly trained within the FAS multigrid framework, and the total loss is PDE residual + load constraint + damping regularization.

Images