Multi-physics field coupling simulation method and device for industrial digital twinning

Abstract

The invention relates to a multi-physics field coupling simulation method and equipment for industrial digital twinning. The method comprises the following steps: automatically dividing a plurality of heterogeneous physical sub-domains based on actual physical field data of industrial equipment, and defining a solving region and a boundary condition for each sub-domain; and then, a variable-driven region division mechanism is utilized to construct a parallel computing map of multiple physical domains, and efficient synchronization of inter-region coupling state variables and real-time interaction of boundary data are realized. In the parallel simulation process, calculation tasks and resources are dynamically and optimally distributed through a self-adaptive scheduling algorithm, bottleneck areas are effectively recognized and eliminated, and the overall simulation efficiency and stability are improved. According to the invention, co-evolution of a plurality of physical processes such as structure, heat, flow, magnetism, electricity and the like can be realized, and unified, flexible and extensible simulation support is provided for digital twinning, state monitoring, performance prediction, control optimization and other applications of a complex industrial system.

Description

Technical Field

[0001] The present invention relates to the field of industrial system modeling and simulation, and in particular to a multi-physics field coupling simulation method for industrial digital twins. Background Art

[0002] Typical industrial scenarios, such as complex equipment manufacturing, power chip thermal control, and intelligent manufacturing control, often involve the state evolution of multiple, cross-coupled heterogeneous physical domains. These include the elastic-thermal response of solid structures, current heat dissipation in power devices, unsteady eddy currents in electromagnetic excitation systems, and heat-fluid transport within coolant channels. Traditional simulation modeling systems typically focus on a single physical field (e.g., FEM structure, FVM heat conduction), making it difficult to achieve a unified expression and coupled solution for state variables across multiple domains at the modeling level.

[0003] Existing multi-domain simulation systems have significant deficiencies in the consistency of data interfaces between models and the uniformity of coupling boundaries. Because multi-domain models often come from different sources and the boundary definition method lacks standardization, logical incompleteness often occurs during region mapping and variable synchronization. Furthermore, traditional task scheduling structures are mostly static and fixed, making it difficult to adapt to the complex requirements of hotspot mutations or bottleneck path transmission during multi-field state evolution. The communication structure between nodes lacks a dynamic reconstruction mechanism, which increases system load and convergence difficulty.

[0004] In multi-domain parallel simulation, it is often impossible to accurately capture the structural bottlenecks and edge weight degradation within the region. Problems such as the attenuation of state information between nodes, imbalance of edge weights, and excessively deep graph structures will directly lead to uneven distribution of computing resources and communication delays. At the same time, many platforms lack a multi-domain linkage scheduling mechanism driven by a graph structure, and are unable to flexibly describe the state flow and boundary transfer between multiple domains in a graph structure, limiting the system's scalability and flexibility. More prominently, there is a lack of efficient fusion scheduling and state monitoring feedback mechanisms between heterogeneous physical domain tasks, making it difficult to achieve system-level task collaboration and regional-level closed-loop solutions, thereby restricting the overall performance and engineering application value of complex multi-physical domain systems.

[0005] Although existing literature has proposed methods such as "graph partitioning scheduling" and "heterogeneous parallel solutions" for multi-task distribution, most of them remain at the task-level load balancing or static grid division, and cannot meet the key demands of industrial digital twin systems in terms of unified expression of cross-physical domain states, information bottleneck identification, graph-level edge reconstruction, and structured scheduling strategy generation.

[0006] Therefore, there is an urgent need for a structural-level coupling modeling and scheduling method that can be used for industrial multi-physical domain systems, support variable sharing between different subdomains, inter-regional state fusion, and system-level communication structure optimization, and provide graph-structure-driven global simulation support for industrial digital twins. Summary of the Invention

[0007] This disclosure provides a multi-physics coupled simulation method and device for industrial digital twins. This method automatically collects and structures multi-source physical field data from actual industrial equipment or systems for the precise construction of coupled simulation models. By intelligently identifying interfaces and coupling boundary conditions between physical domains, it achieves synergy between data-driven and theoretical model-driven approaches.

[0008] One embodiment solves all or some of the shortcomings of model partitioning, regional collaboration, and full-process accelerated simulation in multi-physics industrial systems.

[0009] One embodiment solves all or some of the shortcomings of multi-level scheduling and efficient mapping of computing tasks to achieve the co-evolution of multi-physical domain simulation.

[0010] One embodiment addresses all or some of the shortcomings of efficient boundary reuse driven by cache coherence and prediction at edge nodes in multi-domain coupled simulations.

[0011] One embodiment solves all or some of the shortcomings of graph-based simulation resource scheduling optimization and automatic identification of information transfer bottlenecks.

[0012] One embodiment provides a multi-physics field coupling simulation method for industrial digital twins, including the following steps:

[0013] Constructing a coupled simulation model based on the actual physical field data of industrial equipment or systems, dividing the different physical domains in the coupled model into sub-models, and establishing the solution area and boundary conditions for each physical domain;

[0014] The coupled model is divided into regions using a domain decomposition method to construct a corresponding region-level parallel computing graph, where each region serves as a computing node and the shared boundaries between nodes serve as connecting edges, forming a topological structure for synchronous solution of multiple physical fields;

[0015] Based on the parallel graph, multi-level scheduling and computing task mapping are performed;

[0016] During the simulation process, an edge node cache mechanism is used to store key state variables of each sub-model, and the cached variables are reused in the cross-domain solution synchronization phase to reduce repeated calculations and communication delays;

[0017] Based on the real-time status monitoring results of simulation tasks, an adaptive scheduling diagram is constructed to dynamically optimize simulation resource allocation and solution strategies to achieve accelerated operation of multi-physics field modeling and simulation in industrial digital twin systems.

[0018] According to one embodiment, dividing the coupled model into regions to form a topological structure for simultaneous solution of multiple physical fields includes:

[0019] Construct a variable dependency graph, where vertices represent state variables or physical subdomains, edges represent coupling driving relationships between variables, and edge weights are calculated based on the partial derivatives of the residuals of the control equations with respect to the coupled variables.

[0020] The coupling centrality of each variable is calculated based on the edge weights in the graph, reflecting the strength of its association with other variables;

[0021] Determine the priority of computing nodes based on the coupling centrality results and assign high-coupling nodes to areas with stronger computing power or higher communication bandwidth;

[0022] Establish a scheduling objective function to minimize the cross-region communication cost between highly coupled nodes and optimize the parallel computing topology.

[0023] According to one embodiment, performing multi-level scheduling and computing task mapping includes:

[0024] At the first level, a load balancing strategy is used to build a computation priority scheduling queue based on the computational complexity and communication dependency of the physical fields in each region;

[0025] At the second level, a graph partitioning algorithm is used to map regional computing tasks to processing units or GPU clusters on a multi-core computing platform to form a computing resource assignment table.

[0026] At the third level, based on the time stepping and synchronous triggering mechanism, the iterative process of each area is integrated with tasks, compressed with communications, and scheduled with synchronous solutions to achieve collaborative evolutionary simulation of multiple physical domains.

[0027] According to one embodiment, the first level includes:

[0028] Obtain the coupling centrality and communication complexity indicators of each regional node;

[0029] Calculate the comprehensive scheduling weight of each node to balance computing load and communication intensity;

[0030] Sort all regional nodes according to their weight values ​​and build a scheduling priority queue;

[0031] Each computing task node in the scheduling queue is submitted to the task scheduler in sequence and waits to be allocated to the computing resource unit according to priority.

[0032] According to one embodiment, the third level includes:

[0033] Execute the solution of local coupled control equations in each physical area and complete single-step iterative calculation;

[0034] Compress and transmit boundary variables between regions to reduce the amount of communication data;

[0035] Determine whether the inter-regional coupling residuals meet the synchronous advancement conditions and control the time step consistency;

[0036] If convergence has not occurred, a new round of iterative scheduling is triggered until the co-evolution requirements are met.

[0037] According to one embodiment, the edge node caching mechanism further includes constructing a coupled boundary prediction consistency graph and performing graph-driven cache reuse scheduling, specifically including:

[0038] In each simulation iteration cycle, a cross-region prediction consistency graph structure is constructed based on the coupling boundary relationship between sub-models. Each node in the graph corresponds to a regional boundary cache, and the edge weight represents the correlation of the deviation change trend between the prediction values ​​of two boundary caches.

[0039] For each node, based on its historical cache evolution sequence and the synchronization residual function between adjacent nodes, it iteratively updates its consistency state flag and performs a local consistency consensus voting process to obtain a global decision on whether the node is in the predicted convergence state;

[0040] For all nodes in the consistency graph that are in a converged state, their boundary variables do not require communication synchronization in the next time step, and are directly self-predicted by the edge cache to participate in the simulation solution; otherwise, only the unconverged connected subgraph area in the consistency graph is activated to trigger boundary synchronization and communication update.

[0041] According to one embodiment, the steps of constructing and optimizing the adaptive scheduling graph include:

[0042] Construct a graph structure that integrates physical states and coupling indicators for scheduling optimization;

[0043] Perform multiple rounds of information propagation through graph neural networks to extract key path structural features;

[0044] Identify information transfer bottleneck areas based on information contraction criteria;

[0045] Perform structural fitting analysis and graph reconstruction on bottleneck areas to optimize node connectivity and resource scheduling efficiency.

[0046] According to one embodiment, identifying the information transmission bottleneck area based on the information contraction criterion includes:

[0047] Based on the graph neural network model, multi-level neighborhood information aggregation operations are performed on the adaptive scheduling graph;

[0048] Gradually expand the node structure perception range to obtain the state expression within the multi-hop neighborhood;

[0049] Based on the information contraction theory, the information loss rate of key nodes is calculated to determine whether there is a bottleneck area in the information propagation channel.

[0050] According to one embodiment, performing structure fitting analysis and atlas reconstruction on the bottleneck region includes:

[0051] Select or fit a regular graph structure for the bottleneck area obtained by screening, extract its adjacency matrix and spectral features; calculate the local spectral gap index of the regular graph to measure the tightness of node connections and structural balance;

[0052] When the local spectral gap index is significantly lower than the average level of the entire graph, it is judged that the nodes in this area are not tightly connected and subsequent edge structure optimization processing is required. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] The present invention will be described in more detail below based on embodiments and with reference to the accompanying drawings, wherein:

[0054] Figure 1 The finite element deformation simulation diagram of the rubber material at different loading stages is shown, demonstrating the evolution process of structural stress and boundary deformation.

[0055] Figure 2 Finite volume simulation diagram of current density in the inductor coil area.

[0056] Figure 3 This is the ohmic heat power density distribution diagram in the power device area.

[0057] Figure 4 It is the stress distribution diagram and displacement cloud diagram in the cylindrical structure.

[0058] Figure 5 The simulation results of the heat exchange plate structure and its steady-state temperature field.

[0059] Figure 6 for Figure 5 Finite element simulation results of the steady-state temperature distribution of the heat exchanger plate structure shown.

[0060] Figure 7 A flowchart of a multi-physics field coupling simulation method 100 for industrial digital twins provided by an embodiment of the present disclosure is shown.

[0061] Figure 8A flow chart schematically illustrating a method 300 according to one exemplary implementation of the present disclosure is shown.

[0062] Figure 9 A flow chart schematically illustrating a method 400 according to one exemplary implementation of the present disclosure is shown.

[0063] Figure 10 A typical six-node mesh interconnection structure of computing nodes in some embodiments of the present disclosure is shown.

[0064] Figure 11 A two-dimensional array formed by computing nodes and a multi-hop adjustment path between computing nodes in some embodiments of the present disclosure are shown.

[0065] Figure 12 FIG2 shows a three-dimensional reconstructible connection structure of computing nodes in some embodiments of the present disclosure.

[0066] Figure 13 The figure schematically shows a structural diagram of a data transmission path between computing nodes according to an exemplary implementation of the present disclosure.

[0067] Figure 14 A flowchart of a method 500 according to an exemplary implementation of the present disclosure is schematically shown.

[0068] Figure 15 A flow chart of a method 600 according to an exemplary implementation of the present disclosure is schematically shown.

[0069] Figure 16 The communication connections and contents of edge node caching for different edge nodes using the local priority strategy, the method 600 of the present disclosure, and the centralized strategy are shown.

[0070] Figure 17 A flowchart of a method 700 according to an exemplary implementation of the present disclosure is schematically shown.

[0071] Figure 18 A flow chart of a method 800 according to an exemplary implementation of the present disclosure is schematically shown.

[0072] Figure 19 The multi-hop neighborhood aggregation structure of information propagation and the bottleneck area obtained by screening in some embodiments of the present disclosure are schematically shown.

[0073] Figure 20 A flowchart of a method 900 according to an exemplary implementation of the present disclosure is schematically shown.

[0074] Figure 21The diagram shows a ring-of cliques regular graph and an equidegree random graph generated by fitting an original graph that does not meet the regular graph standard using the ring-of cliques method and the equidegree random method in some embodiments of the present disclosure.

[0075] Figure 22 The node arrangement of rows and columns in the adjacency matrix of the original graph provided in an embodiment of the present disclosure is shown.

[0076] Figure 23 FIG. 1 is a block diagram of an electronic device 10 according to an embodiment of the present disclosure. DETAILED DESCRIPTION

[0077] The following description of exemplary implementations of the present disclosure is provided in conjunction with the accompanying drawings, which include various details of the exemplary implementations of the present disclosure to facilitate understanding. These details should be considered merely exemplary. Therefore, those skilled in the art will recognize that various changes and modifications may be made to the exemplary implementations described herein without departing from the scope and spirit of the present disclosure. Similarly, for the sake of clarity and conciseness, descriptions of well-known functions and structures are omitted in the following description.

[0078] As used herein, the term "including" and its variations represent open inclusion, i.e., "including but not limited to." Unless otherwise stated, the term "or" means "and / or." The term "based on" means "based at least in part on." The terms "an example exemplary implementation" and "an exemplary implementation" mean "at least one example exemplary implementation." The term "another exemplary implementation" means "at least one other exemplary implementation." The terms "first," "second," etc. may refer to different or the same objects. Other explicit and implicit definitions may also be included below.

[0079] With the rapid growth in industrial system complexity and multi-field coupling, traditional simulation modeling methods based on a single physical domain or control strategy are no longer able to meet the real-time and precision requirements of digital twin systems in structural prediction, performance evaluation, thermal regulation, electromagnetic response, and energy optimization. For example, high-end manufacturing equipment, energy systems, intelligent machine tools, power electronics modules, and precision structural components are no longer systems dominated by a single physical process, but rather highly complex, cross-domain coupled systems with the following characteristics:

[0080] Ⅰ) There is an interaction between solid structure and thermal expansion (such as Figure 1 rubber materials in the

[0081] II) There are thermal-electro-magnetic response behaviors of power devices (such as Figure 2 Coil, Figure 3 power components);

[0082] III) There is energy and momentum exchange between the structure and the cooling fluid (e.g. Figure 5 – Figure 6 shown);

[0083] IV) There is coupled migration of stress, electric field and heat flow between materials (such as Figure 4 structural response);

[0084] These physical processes occur at different locations, at different time sequences, and at different scales, but must be considered holistically in industrial systems. Otherwise, misjudgments, biased interpretations, and even equipment failure prediction failures will occur.

[0085] The goal of this invention is not simply to "group models from different application fields together", but to build a unified simulation framework across multiple engineering subsystems with real industrial systems as the target object, to realize digital twins with the ability to co-evolve multiple physical processes such as "structure-field-flow-electromagnetism-magnetism".

[0086] For example:

[0087] A. A power module cooling system includes both power chip current distribution (such as Figure 3 As shown), the coil excitation electron cloud distribution (as shown Figure 2 As shown), the temperature gradient of the heat exchange plate (as shown Figure 5 ), and also includes a fluid heat exchange path (as shown Figure 6 shown);

[0088] B. A rubber mechanical trigger or electromechanical actuator under high heat and high pressure ( Figure 1 ) transmits the response signal to the electromagnetic excitation system ( Figure 2 ) or heat transfer structures (such as Figure 5 shown);

[0089] C. Structural parts in industrial systems (such as Figure 4 As shown in the figure, the material may be in a coupled state of current field, temperature field and stress field at the same time, such as in the process of hot rolling, welding and induction hardening.

[0090] Although these scenarios belong to "seemingly unrelated" fields, what they have in common in industrial systems is that the physical processes are "truly" coupled and cannot be decoupled.

[0091] Therefore, only by graph-structuring, unified modeling, and parallel simulation of these essentially interrelated physical processes can we truly build a "full life cycle simulation mapping engine" for industrial digital twins.

[0092] Figure 7A flowchart of a multi-physics field coupling simulation method 100 for industrial digital twins provided by the present disclosure is shown, comprising the following steps:

[0093] 1) Build a coupled simulation model based on the actual physical field data of industrial equipment or systems, divide the different physical domains in the coupled model into sub-models, and establish the solution area and boundary conditions for each physical domain;

[0094] In step S102, including S202: adopting strongly coupled finite volume-finite element hybrid modeling of multi-physics coupled system:

[0095]

[0096] Where u is the displaced mass and γ is the stress tensor; Represents the gradient due to temperature T The thermal expansion stress term caused by elec (E ele ) is due to the electric field strength E ele The electric potential field caused by ele =-φ, Represents the gradient calculation of the corresponding physical quantity, Q joule is Joule heat, Q joule =σ|E ele | 2 , Q mech is the internal heat generated by material deformation; c p is the specific heat capacity at constant pressure, with the unit of J / (kg·K). It determines the speed of material temperature rise and affects the temperature-time dynamic response in heat conduction and thermal stress simulation. ∈ is the dielectric constant in the electric field Poisson equation, with the unit of F / m (farad / meter). It is used to measure the material's response to the electric field and determines the electric field penetration. ρ e is the free charge density (i.e. the net charge per unit volume), with the unit of C / m 3 (Coulomb / cubic meter), which exists as an electric field source term and determines the bias of the electric field distribution; it is the main variable in scenarios such as capacitor charging, semiconductor excitation, and plasma discharge. Calculates the gradient function.

[0097] S204: Hybrid region modeling and partitioning. Within the simulation domain Ω, define subdomains:

[0098] For the simulation domain Ω, N discrete physical domains are used for boundary variable transfer, and the Nth physical domain is Ω (N) , And satisfy Ω (i) ∩Ω (j) =Γ ij , Γ ij is the i-th physical domain Ω(i) The boundary between the jth physical domain; i = 1, 2, ..., N; it is only necessary to set the boundary condition type (such as temperature, normal heat flux, electric potential, current density, displacement or stress) of the corresponding regional mesh during the regional mesh generation stage.

[0099] If the FEM+FVM hybrid method is used for discretization, define Ω (1) To use the finite element method (FEM) for the elastic and thermal stress regions, define Ω (2) Finite volume model (FVM) is used for strong heat transfer and electric potential distribution area (also called strong thermal conductivity and electric potential distribution area); the subdomains are separated by interfaces. Perform border exchange.

[0100] By constructing the local discretized equations as follows:

[0101] K (1) u (1) =F (1) ;

[0102]

[0103] G (2) φ (2) =b (2) ;

[0104] The first equation in this discrete set of equations is the modeling equation for the elastic and thermal stress domain (FEM), the second equation is the modeling equation for the thermal conductivity domain (FVM), and the third equation is the modeling equation for the electric potential domain (FVM). The brackets and numbers in the upper right corner of this discrete set of equations represent the simulation domains where the modeling equations formed by different methods are located. The (1) in the upper right corner indicates that the variable / matrix belongs to the simulation domain Ω. (1) , the domain uses the finite element (FEM) method; the upper right corner is (2) indicating that the variable / matrix belongs to the simulation Ω (2) , the finite volume method (FVM) is used in this domain. The first equation states: After discretization by FEM, in the structural domain Ω (1) Internal, nodal displacement u of the elastic structure (1) Under force F (1) Under the action of stiffness matrix K (1) determined;

[0105] The second equation states: Heat conduction in Ω (2) The time-varying dynamic response of the heat capacity C (2) and thermal conductivity K (2) control; is the temperature change rate.

[0106] The third equation expresses the potential distribution in the electric field governed by the Poisson equation, represented by the conductivity matrix G. (2) and source term b(2) Decide.

[0107] The remaining discrete physical domains can be electromagnetic coupling domains (such as induction heating, EM field), and the second physical domain - electric field Ω (2) and the first physical domain Ω (1) Forming thermal-electromagnetic coupling; it can also be a flow field physical domain (such as coolant channel, natural convection), and Ω (1) Forming fluid-solid coupling (wall shear and heat transfer), and at the same time with Ω (2) Form a temperature field connection; it can also be a chemical reaction subdomain (such as electrode, fuel cell, hydrogen evolution, etc.), which can be connected with Ω (2) The electric field and temperature field form thermal-electro-chemical coupling.

[0108] The above model ensures that each subdomain uses the optimal method to construct the numerical structure, thereby improving the solution efficiency.

[0109] Figure 1 The finite element simulation diagram of the rubber material at different loading stages shows the elastic structure subdomain Ω (1) The local deformation and stress-strain evolution responses during the initial compression, intermediate loading, and final loading processes intuitively reflect the expansion trend of the strain area along the heat source direction and its relationship with the surrounding subdomains (such as the fluid domain Ω (5) ) between the boundary shape variables. Figure 1 (a)– Figure 1 (b) shows Ω (1) Typical simulation results under the domain (elastic structure + thermal expansion stress field), specifically Figure 1 (a) shows the elastic structure subdomain Ω (1) Local deformation response at the initial stage of compression. Figure 1 (b) shows the sub-domain Ω (1) The process of the medium strain region expanding along the direction of the heat source. Figure 1 (c) shows the Ω (1) and surrounding subdomains (such as fluid domain Ω (5) ) between the boundary shape matching behavior. Figure 1 The boundary variables of the middle domain include displacement u, stress tensor γ, normal force and heat flux input. ( is the temperature gradient coefficient), which is directly dependent on the temperature field gradient. Therefore, if the deformation distribution of the structural area in the figure is obviously affected by the hot zone, it can be judged that the deformation distribution is related to Ω. (2) There is a shared temperature boundary in the domain; if the deformation region occurs near the fluid cooling surface, it can also be determined that (5) There is a wall heat flux boundary in the domain, and the graph structure is represented by Ω (2) Domain, Ω (5) Domain pointing to Ω (1) Create directed edges across domains.

[0110] Figure 2 、 Figure 5 and Figure 6 Ω (2) Finite volume analysis diagram related to the domain (region of strong thermal conductivity and electric potential distribution), with boundary conditions including temperature T, heat flux q, electric potential φ, and current density J. Figure 2 Schematic diagram of the electric field subdomain Ω (3) In the field of electric heat conduction Ω (2) Current accumulation paths formed by the interaction of boundary conditions. Figure 5 For the simulation result of the steady-state temperature field, if the thermal field and electric field (such as Figure 2 Ω (3) domain) or fluid region (such as Figure 6 The fluid cooling domain, Ω (5) domain) has boundary continuity or transfer characteristics, it can be judged that there are multiple outward directed edges on it; specifically, Ω (2) As an energy source or potential coupling intermediary, it plays the role of "source node" in multiple physical domains in the graph structure. Figure 5 (a) shows a heat exchange plate model for heat dissipation of electronic equipment. The inlet cooling medium is pure water, the volume flow rate of the cooling medium at the inlet is set to 8L / min, and the inlet cross-sectional area is 80mm 2 , the corresponding average flow rate of the cooling medium at the inlet is 1.67m / s. The heat exchange plate material is 6063 aluminum alloy, and the flow channel surface area A=60,000mm 2 , set the heat load P = 1500W, the inlet cooling medium temperature is 18 ° C. The outlet temperature T out The state of the fluid in the flow channel is automatically solved by the simulation process.

[0111] Figure 5 (b) shows Figure 5 (a) The finite element simulation results of the steady-state temperature distribution of the heat exchange plate show the heat conduction subdomain Ω (2) With the fluid area Ω (5) Heat is exchanged between them through the boundary temperature difference. Figure 6 (a) shows the coolant velocity field distribution diagram in the cooling medium channel of the heat exchange plate. Figure 6 (b) shows the coolant temperature field distribution diagram, which is used to characterize the fluid field Ω (5) With solid domain Ω (1) and thermal domain Ω (2) Multi-boundary coupling behavior between.

[0112] Figure 2 The current contour diagram results of a high current density area concentrated in a coil (a coil is also a power component) are also shown, and its boundaries include the voltage difference between conductors, charge source and equipotential surface settings.

[0113] Figure 3 is the ohmic heat power density distribution diagram in the power device area, Figure 3 The figure shows the internal ohmic heat distribution of an induction coil winding (which belongs to a power device), namely Figure 3 shows the electromagnetic induction domain, Ω (4) Domain, schematically showing the electromagnetic induction subdomain Ω (4) Heat source distribution and its relationship with electric field Ω (3) The boundary source coupling relationship between them. Figure 3 The ohmic heat distribution inside the induction coil winding is Figure 2 The current density in the space is dual. If the ohmic heat is mainly affected by Ω (3) Domain or Ω (2) The electric field control effect of the domain, then Ω (4) The domain needs to be from Ω (2) The frequency excitation current source data is obtained from the domain and an input edge is established in the graph structure. Figure 3 The structure is located around the heating area and can also determine Ω (4) Domain and Ω (1) There is a magnetic-thermal boundary in the domain.

[0114] Figure 6 A fluid velocity and temperature field distribution is shown. The flow velocity and heat-flux distribution are solved by FVM. The boundaries include inlet flow velocity, wall shear stress, and convective heat transfer coefficient. Figure 4 It further shows the effect of coolant heat exchange on the cylindrical structure. Figure 4 is the stress distribution diagram in the cylindrical structure ( Figure 4 (a)) and displacement cloud map ( Figure 4 (b)) is used to illustrate the structural field Ω (1) In the thermal basin Ω (5) and the potential domain Ω (2) The composite response under double boundary action; thus, it can be confirmed that Ω (5) Respectively with Ω (1) (fluid-solid boundary), Ω (2) (Temperature source boundary) There are at least two boundary interactions. In the graph structure, Ω (2) and Ω (1) Unidirectional Ω (5) Sending heat flux or temperature gradient information can construct two points in the graph structure to Ω (5) The directed edge of .

[0115] The above five subdomains pass through the interface boundary Γ ij =Ω (i) ∩Ω (j)Perform row-by-row interaction coupling and mark the required boundary exchange variables, such as potential-heat, heat-flow, stress-displacement, etc., during the initial mesh generation stage to ensure boundary consistency and solution closure under hybrid modeling.

[0116] As another embodiment of the present invention, Figure 8 A flowchart of method 300 is shown schematically illustrating an exemplary implementation of the present disclosure, which uses a domain decomposition method to divide the coupling model into regions and construct a corresponding regional parallel computing graph, wherein each region serves as a computing node and the shared boundaries between nodes serve as connecting edges to form a topological structure for synchronous solution of multiple physical fields.

[0117] In the region division and parallel graph construction steps, the simulation area will be topologically divided according to the "coupling variable influence path" instead of simply dividing the area according to the physical location.

[0118] In step S302, a variable-driven dependency graph G is constructed. v =(V,Y)

[0119] Vertex v i Represents the i-th state variable or physical subdomain (such as temperature field unit, electric potential field unit);

[0120] Edge Y ij Represents a coupling driving relationship between variables i→j (such as temperature gradient Influence stress tensor γ);

[0121] The edge weight can be set as a sensitivity index of the variable coupling coefficient: where R j is the residual of the j-th variable control equation, x i is the i-th coupling variable. The larger the weight, the stronger the interdependence, and the more important it is to arrange them in parallel.

[0122] In step S304, the parallel partitioning scheduling strategy is adjusted to "coupling centrality driven": the coupling centrality C of each subdomain or each variable is calculated. i , as a priority indicator for parallel computing:

[0123]

[0124] In step S306, high C i Nodes are preferentially allocated to core nodes or high-bandwidth computing resources to ensure minimum delay in the coupling transmission path.

[0125] In step S308, the scheduling objective function is:

[0126]

[0127] That is, nodes v that communicate frequently are i , v j Mapped to the same processor, reducing cross-unit communication costs (such as PCIe transmission or inter-cluster communication). ij and represent the coupling influence weights between node i and node j, respectively, which measure x i R j sensitivity; w ji Represents the coupling influence weight between node j and node i, which measures x j R i sensitivity; w ji represents the coupling influence weight between node j and node i; M(v i )≠M(v j ) indicates the communication cost M(v i ) and the communication cost M(v j ) are not equal;

[0128] This minimizes the weighted communication cost after coupling, ensuring that highly coupled variables are distributed to topologically adjacent computational units whenever possible. CommCost(i,j) represents the communication cost between nodes i and j. By constructing a dependency graph based on the influence relationships between coupled variables, the multiphysics model is partitioned into regions based on the coupling strength rather than physical location. Compute node priorities are determined based on variable coupling centrality, optimizing the parallel computing scheduling structure to reduce cross-node communication overhead.

[0129] Complete the above Figure 1-5 、 Figure 2-Figure 3 、 Figure 4 as well as Figure 5-Figure 6 After the physical domains are partitioned, the system enters the parallel scheduling structure construction phase. According to steps S302–S308, the five subdomains are mapped in parallel using a variable-driven domain partitioning method to construct a regional computational graph, where each subdomain corresponds to a computational node (or subgraph). The edges in the graph represent the boundary interactions between domains.

[0130] Figures 1 to 6 They respectively reflect the computational complexity of different physical domains, the boundary exchange strength and the coupling relationship in the graph structure. For example, Figure 2 Shows the distribution of concentrated areas of high-intensity current density, representing the electric field Ω (3) The solution result is Depends on the electric potential field f elec (E ele ) distribution gradient, so it is necessary to start from the electrothermal coupling area Ω (2) Receives boundary potential input. Figure 2 The region shown also contains Ω(2) Domain and Ω (3) domain, both in the boundary region Γ 23 When the corresponding regional parallel computing graph is constructed in step S302, the computing graph structure constructed is composed of Ω (2) Domain pointing to Ω (3) Domains establish directed edges, indicating that variables are transferred from the former to the latter.

[0131] Another example Figure 6 The upper and lower figures respectively represent the distribution of velocity field and temperature field in the coolant channel. Figure 6 In the figure above, the velocity gradient mainly appears in the bend and transition areas of the flow channel, reflecting the fluid area Ω (5) physical behavior; and Figure 6 In the figure below, the coolant temperature shows a distribution gradient that clearly matches the heat source area, indicating that the fluid domain receives heat from Ω on the boundary. (2) The temperature input of the (heat conduction domain) shows a typical thermal-fluid boundary coupling phenomenon. Therefore, when performing step S302, the Ω in the constructed calculation spectrum structure is (5) Domain and Ω (2) There is a boundary Γ between domains 25 Therefore, in the graph structure, Ω (2) Domain and Ω (5) A directed edge representing heat transfer should be established between domains.

[0132] The multi-physics field coupling simulation method 300 for industrial digital twins provided in this embodiment can effectively reduce the communication distance between highly coupled variables, reduce global synchronization consumption, and improve the efficiency and scalability of system parallel simulation. It is particularly suitable for large-scale, strongly coupled multi-domain coupling computing tasks.

[0133] As another embodiment of the present invention, Figure 9 The flowchart of a method 400 for performing multi-level scheduling and computing task mapping based on the parallel graph according to an exemplary implementation of the present disclosure is schematically shown, specifically comprising:

[0134] a. Figure 9 The specific steps of building a computation priority scheduling queue using a load balancing strategy at the first level are shown according to the computational complexity and communication dependency of the physical fields in each region.

[0135] Specifically, in step S402, during the scheduling preparation phase, the coupling centrality C of each computational region node must be extracted. i and communication complexity Comm i Two indicators; among them, coupling centrality C iIt is used to measure the physical interaction strength between the node and other coupled subdomains, and has been obtained based on the edge weight calculation in the variable dependency graph; the communication complexity Comm i It reflects the data exchange frequency and communication load of the node in the parallel graph, which can be obtained by counting the shared boundary length between it and the adjacent nodes, the data exchange frequency, the historical communication overhead or the amount of data interaction per unit time; the above two indicators correspond to the computing density dimension and the communication sensitivity dimension respectively, and jointly constitute the basis of the subsequent scheduling priority scoring function.

[0136] In step S404, combined with the coupling centrality C i , communication complexity Comm i , and calculate the comprehensive scheduling weight: W i =αC i +βComm i , where α+β=1; D i It is usually determined by factors such as the number of internal grids, the order of the control equations, and the number of time step iterations; Comm i represents the communication dependency index between the i-th computational region and other regions, which can be represented by the shared boundary length, communication frequency, or data interaction volume. α and β are artificially set weight parameters used to balance the computational workload and communication workload (e.g., computation-intensive: α = 0.7; communication-sensitive: β = 0.6).

[0137] Communication complexity Comm i By quantifying each regional node v i It is defined by the boundary communication load of all its adjacent regional nodes and is expressed as:

[0138] in, For the i-th node v i The index set of all adjacent nodes with boundary interactions; L ij For the i-th node v i With the jth node v j The length of the shared boundary or the number of boundary grids between them; μ ij is the data exchange frequency or communication load factor per unit shared boundary per unit time, which can be determined based on simulation history statistics or the average size of data exchanged per step; if time correlation and compression effects are considered, then is the average amount of data transmitted per unit time, ρ ij ∈(0,1] is the compression ratio factor (1 if not compressed).

[0139] Communication complexity Comm iThe larger the value, the higher the data interaction volume of the node in the multi-physical domain coupling. When scheduling, it should be allocated to a location with sufficient bandwidth resources or closer to high-frequency interaction nodes. If the graph is an unweighted graph, it can be simplified to count the number of adjacent edges and the average interaction data intensity. In parallel scheduling, Comm i It is often used as a communication cost evaluation factor in the task mapping objective function to avoid communication bottlenecks.

[0140] In step S406, according to the weight W i Prioritize the nodes to form a scheduling priority queue (hereinafter referred to as the scheduling queue):

[0141] {v i} represents the set of all regional nodes, where v i is the i-th computational sub-area node in the graph, i.e. the i-th physical domain Ω i For the corresponding computing task node, the i-th physical domain is the i-th computing sub-area; To schedule a computing task priority queue, press W i Arrange the node set from high to low {v i}; Each They will enter the scheduler in order of ranking and wait for allocation to specific hardware resources.

[0142] b. At the second level, a graph partitioning algorithm is used to map regional computing tasks to processing units or GPU clusters on a multi-core computing platform, forming a computing resource assignment table;

[0143] According to the graph partitioning method, graph J is divided into K subgraphs:

[0144]

[0145] Each subgraph is mapped to a processor core or GPU cluster P k , generate computing resource allocation table:

[0146] Figure J is the original task graph, including all regional nodes v i and the communication edges between them; J i Including some nodes and edges, K represents the total number of processor cores or GPU clusters, that is, the number of parallel hardware resources; P k The kth processing unit (such as CPU core / GPU block) is used to execute the tasks in a subgraph; Load(J i ) represents the total computational load in the i-th subgraph, such as the number of grids, the number of task steps, or the estimated floating-point operations FLOPs; Balance means Load(J1)≈Load(J2)≈…≈Load(JK );

[0147] The goal of this condition is to ensure that each processing unit is assigned a roughly equal amount of work, avoiding load imbalance that causes idle resources or task bottlenecks.

[0148] M(v i ):v i →P k , which minimizes the cross-node communication cost. This formula is the resource mapping function, which means that the task node v i Assigned to hardware unit P k 。 M(v i ) represents the i-th node v i communication cost.

[0149] As the computational complexity of multi-domain coupled simulation continues to increase, a single computing node can no longer meet the timeliness requirements of high-precision field coupling models under large-scale parallel conditions. Therefore, the parallel scheduling method using regional division and computing node graph organization structure has gradually become the mainstream solution. In order to adapt to the differences in coupling strength and communication structure requirements between different types of physical domains, Figure 10 、 Figure 11 and Figure 12 Different scheduling and connection methods of computing nodes in the present disclosure are schematically shown respectively. Figures 10 to 12 The following diagram schematically shows a typical scheduling structure diagram of multiple computing nodes under different connection topologies, where each node can contain multiple computing devices, and the nodes can interact with each other through different boundary connection relationships. For example, Figure 10 shows a typical six-node mesh interconnection structure, Figure 11 shows a multi-hop scheduling path between computing nodes in a two-dimensional array, Figure 12 A three-dimensional reconstructed connection structure is shown. These topological structures can serve as the basis and reference examples for constructing data transmission paths between nodes in task scheduling graphs.

[0150] Figure 10 The figure shows the central cross topology connection structure between the six queues, where nodes A, B, C, D, E, and F represent the six computational sub-areas constructed in steps S302–S308, and each area is decomposed by the variable coupling path. Figure 16 Node C, as the node with the highest coupling centrality, has strong coupling relationships with the other four nodes, forming a cross-connected structure. Each node is shown below its contained computational unit groups (such as T1–T8, U1–U8, V1–V8, W1–W8, H1–H8, and R1–R8). These computational units are used to perform the discrete solution of the local governing equations in parallel.

[0151] Figure 11This diagram schematically illustrates a two-dimensional grid connection structure suitable for simulation scenarios where the coupling between sub-regions is relatively balanced and the computational tasks are regularly distributed. In this structure, nodes are arranged in a two-dimensional array and connected by adjacent edges, forming a task partitioning map driven by spatial topology. This structure is particularly suitable for modeling and solving problems where the physical domain itself has a regular grid structure (such as two-dimensional heat conduction and multilayer thin film materials).

[0152] Figure 12 A graph structure combining a deformed cube grid and redundant surrounding edges is presented, suitable for modeling multi-level coupling problems in heterogeneous scenarios. In this structure, five nodes are arranged in three dimensions via multiple redundant connection paths, while the multiple computational sub-units contained within form a locally coupled subgraph. In practical implementation, this structure improves scheduling redundancy, reduces bottleneck edge bandwidth load, and provides structural support for subsequent graph reconstruction and edge flipping optimization steps (such as information bottleneck identification and the edge adjustment mechanism in S706–S708).

[0153] Overall, Figure 10 – Figure 12 This constitutes an illustration of three typical task scheduling topologies supported by the variable-driven region partitioning method disclosed in the present invention (a method for constructing a task scheduling topology structure between sub-regions in a computational graph in steps S302–S308). The appropriate graph structure can be adaptively selected based on the coupling characteristics of a specific simulation model and deployed in a parallel computing system.

[0154] In addition, the present disclosure provides a method for constructing a computing array based on graph structure partitioning and task scheduling optimization, which realizes a task partitioning and scheduling strategy with more topological cognition capabilities by constructing a coupling graph structure between nodes. Figure 13 , further explains the data transmission path relationship between computing nodes.

[0155] Figure 13 The diagram schematically shows a data transmission path structure diagram between computing nodes (processing units or GPU clusters of a multi-core computing platform) according to an exemplary implementation of the present disclosure. The figure shows a computing array of 8 rows and 8 columns, corresponding to 8 computing nodes (nodes Ndoe1 to Node8, i.e., 8 processing units or 8 GPU clusters), each computing node includes 8 computing devices (i.e., each processing unit or each GPU cluster includes 8 computing devices), numbered V11 to V88 respectively. Each row represents a computing node, and the 8 computing devices in each node are arranged side by side in a horizontal direction. Through this structure, a processing unit 200 or a GPU cluster 200 can be formed for performing discrete solutions of multi-physics control equations in parallel.

[0156] According to an exemplary implementation of the present disclosure, the computing devices in the same computing node can come from the same or different physical device deployment environments. For example, V11 to V18 can all be deployed on the same GPU hardware, or they can be distributed in two independent GPU systems. Furthermore, depending on the functional expansion capabilities, the gray diagonal filled area on the right side of the figure indicates scalable computing devices (such as V17, V18, V27, V28, etc.), which can support node expansion when the subsequent simulation scale is expanded, such as by horizontally splicing new node arrays to expand their physical partitions or functional dimensions.

[0157] like Figure 13 As shown, the figure draws a representative data transmission path connection. Each arrow represents the transmission of state variables or boundary control information from a computing device in a certain node to a target device in another node. Taking the path (V11→V44) as an example, it means that there is a high coupling relationship between the device D11 in the node Node1 and the device V44 in the node Node4. In this coupling graph, the starting and ending nodes of the path are determined by sorting the node set with a larger comprehensive scheduling weight in the variable dependency graph in sequence according to the scheduling queue, which belongs to the specific implementation form of the specific configuration of multiple processing units in step S406 according to the computing resource assignment table. The number of paths shown in the figure is a simplified typical connection, which is used to represent the state information transmission process between physical coupling areas.

[0158] It should be noted that the above paths support bidirectional communication, that is, the boundary values ​​of the state variables can be provided by the sending node or fed back by the receiving node, which is used to achieve regional synchronous update during the time step. In addition, the connections between nodes in Figure x follow the principle of minimizing the communication cost of the coupling centrality path, and the path layout follows the weight W proposed in S404-S406. i The scheduling queue formed by prioritizing nodes is to map nodes that communicate frequently to topologically adjacent areas as much as possible, thereby reducing communication delays and network congestion.

[0159] pass Figure 13 The processing unit or GPU cluster structure of the multi-core computing platform for receiving the regional computing tasks mapped by the graph partitioning algorithm shown in the figure only needs to reserve a single port for each computing device to access the cross-node communication, thereby saving port resources in hardware deployment. Figures 10 to 12The internal loop connection is constructed in any of the ways described in , and is not limited to a specific layout, allowing dynamic composition and task migration. Furthermore, the computing array structure has good scalability and adaptive scheduling capabilities, and is suitable for multi-point data synchronization operations in computing scenarios (such as tensor summation, gradient reduction, etc.), which can significantly improve the communication efficiency and load balancing in large-scale parallel simulations. Through the method 400 disclosed in the present invention, it is possible to dynamically adapt to different computing resources and task distributions, improve the balance of task allocation, give full play to the potential of hardware resources, and at the same time reduce the blocking probability of high-load / high-communication nodes, thereby significantly improving the overall system throughput and computing concurrency capabilities.

[0160] Figure 14 A flowchart of method 500 is shown schematically illustrating a third level of an exemplary implementation of the present disclosure for performing task fusion, communication compression, and synchronous solution scheduling on the iterative process of each region to achieve collaborative evolutionary simulation of multiple physical domains based on time stepping and synchronous triggering mechanisms.

[0161] In step S502, for the i-th physical domain Ω i Solve the local strongly coupled governing equations:

[0162]

[0163] Among them, The residual expression of the strongly coupled field control equations constructed in step 1) in the ith physical domain at the t+1th time step, in which the values ​​passed in from the adjacent boundaries in the previous time step are used It is solved as boundary conditions; the solution method can be implicit Newton iteration method or relaxation iteration.

[0164] In step S504, at the boundary Γ between the i-th physical domain and the j-th physical domain, ij Variable compression transmission is performed between nodes i and j, that is, when there is boundary coupling between nodes i and j, the transmission variables on the shared boundary are interpolated and compressed:

[0165]

[0166] is the compressed boundary vector at the t+1th time step (used to reduce the amount of communication data); InterpCompress(·) is the boundary variable interpolation compression function, which can be implemented through polynomial interpolation (such as cubic Lagrangian), wavelet compression (such as Haar), or boundary mapping compression after dimensionality reduction based on the POD mode.

[0167] The boundary vector after interpolation compression Propagate to adjacent physical domains for boundary update in the next iteration.

[0168] In step S506, the convergence of the coupling residuals between regions is checked and the time step is uniformly advanced. Specifically, the global synchronous controller performs a coupling consistency check to determine whether all regions have met the synchronous convergence conditions:

[0169] Define the region Ω at the t+1 time step i and Ω j The shared region boundary Γ ij The coupling residual on :

[0170]

[0171] If all residuals satisfy: Then all computing nodes advance to the next time step t→t+Δt;∈ tol is the convergence residual iteration termination threshold;

[0172] Otherwise, return to step 1) and perform the next round of iteration (using the Jacobi iteration method or the Gauss-Seidel iteration method).

[0173] At the t+1th time step, the region Ω i The calculated state variable u i , in the neighborhood Ω j Bordering Boundary Γ ij This is equivalent to solving the problem inside the region. It is "truncated" or "extracted" at the boundary to participate in boundary matching or compare with neighborhood values.

[0174] 4) During the simulation process, an edge node cache mechanism is used to store key state variables of each sub-model, and the cached variables are reused in the cross-domain solution synchronization phase to reduce repeated calculations and communication delays;

[0175] 5) Based on the real-time status monitoring results of simulation tasks, an adaptive scheduling diagram is constructed to dynamically optimize the simulation resource allocation and solution strategy to achieve accelerated operation of multi-physics field modeling and simulation in the industrial digital twin system.

[0176] Through the method 500 disclosed in the present invention, the bandwidth pressure of high-frequency boundary data transmission in multi-domain simulation is effectively reduced, the rapid convergence of cross-domain synchronization under large-scale parallel computing is achieved, the coordination of solutions between physical domains and the simulation stability are improved, and technical support is provided for high-precision engineering simulation and real-time scenarios.

[0177] Figure 15 A block diagram schematically illustrates a process 600 of constructing a coupling boundary prediction consistency graph and performing graph-driven cache reuse scheduling according to an exemplary implementation of the present disclosure.

[0178] Specifically, in some embodiments, in step S602, within each simulation iteration cycle, the cross-region prediction consistency graph structure constructed based on the coupling boundary relationship between the sub-models is the variable-driven dependency graph G constructed in step S302. v =(V,Y), the i-th node v in the graph i With the jth node v j The boundary Γ ij The cache is a region boundary cache corresponding to each node in the graph, R j is the residual of the j-th variable control equation, then the edge weight w ij Then it is the correlation of the deviation change trend between the two boundary cache prediction values. In addition, define x i ∈{0,1},x i For the i-th node v i Convergence state prediction value (which is a Boolean variable) that indicates the node v i Is it in the prediction convergence state (if x i =1, then node v i is a convergent state; if x i =0, it means that node v i is not converged), define y ij ∈{0,1},y ij is the reference value of the edge residual consistency constraint (which is a binary ratio), if y ij =1, indicating that node v i With adjacent node v j The edge Γ ij The residual consistency carried is turned off, that is, the synchronization residual on this edge has a constraining significance for the decision; if y ij = 0, then edge Γ is not considered ij , that is, the adjacent node v j The information of node v i The state determination of the adjacent node v can be ignored. j The information is invalid or node j has converged.

[0179] Define z i is the historical cache prediction residual (which is a continuous variable), z i =||ΔR i (t)||2, that is is the i-th node v of the t-th generation simulation round i The qth sample value of the residual of the variable control equation, is the i-th node v in the t-1 generation simulation round i The qth sample value of the residual of the variable control equation, d is the total number of samples; q = 1, 2, 3, ..., d;

[0180] In step S604, an objective function is constructed to minimize the total influence of the edge weight residual of the non-converged subgraph in the graph:

[0181] When x i x j =1, indicating that the i-th node v i and the jth node v j If both converge, then y ij The previous coefficient is 0, so you can skip y ij constraints;

[0182] a. Construct local convergence residual control constraints:

[0183]

[0184] Define the i-th node v i Convergence state prediction value x i The piecewise linear function of i (x i ), Θ i (x i )=||ΔR i (t)||2-ε·x i , when x i = 0, Θ i (x i ) is ||ΔR i (t)||2; when x i =1, Θ i (x i ) is ||ΔR i (t)||2-ε;

[0185] The local convergence residual control constraint of formula (A1) indicates that the residual should be lower than the threshold ε when the prediction converges, otherwise it is forced to 0. That is, the local convergence residual control constraint of formula (A1) is used to determine whether node i is in the prediction convergence state. When x i =1, the residual is forced to be less than the threshold ε, otherwise there is no restriction.

[0186] b. Constructing neighborhood consistency coupling constraints:

[0187]

[0188] That is, for any node v that belongs to V i , node v j , both have This constraint condition holds; α ij represents the predicted stability influence of node j on node i, and β is the constant coefficient for suppressing invalid edges.ij and influence coefficient α ij , reflects the impact of whether the adjacent nodes converge on the lower limit of the current node residual, reflecting the cross-node state coordination. If the adjacent node j does not converge (at this time x j =0) and y ij =1, then the constraint is activated, requiring the residual z of this node i It cannot be less than this influence value; otherwise, the constraint is "soft closed".

[0189] c. Construct subgraph voting consensus constraints:

[0190]

[0191] In this constraint, Vote i (y ij ) is the bid support function (i.e., the function for calculating the number of consistent votes obtained). This constraint indicates that if node i wants to be judged as "locally consistent", it needs to obtain at least The voting support of neighboring nodes is used to trigger the consensus. Only when most of the neighboring "identifying" nodes in the neighborhood are in the same state can it be judged as local convergence. Therefore, the subgraph projection consensus constraint expressed by formula (A3) is equivalent to a local majority consensus mechanism. This constraint ensures that each node participates in the voting process in the local neighborhood.

[0192] d. Construct boundary communication activation switching constraints:

[0193]

[0194] That is, for all i-th nodes v i , all satisfy and This constraint is used to activate communication boundary synchronization, indicating whether the communication synchronization path (non-converged subgraph) is activated. G is the activation indicator variable for communication synchronization. When G = 1, it indicates that the current region has not converged and the boundary needs to be synchronized. When G = 0, it indicates that the current region has converged and synchronization is not required.

[0195] For the i-th node v i The corresponding residual consistency constraint coefficient is used to allow the prediction residual to be lower than the threshold, which constrains the formula (A1); For the i-th node v i With the jth node v j The synchronous residual coupling term is used to smooth the activation of the neighborhood constraint, which constrains the formula (A2); is the communication activation control coefficient, which is used to indicate whether the boundary triggers synchronization control, and constrains formula (A4); and They are all slack variable coefficients, and the control constraint boundaries are within the feasible physical domain.

[0196] In step S604, a local consistency consensus voting process is further performed to obtain a global decision on whether the node is in the predicted convergence state. Specifically, the interior point method barrier function target model is constructed:

[0197]

[0198] Among them, the first is the actual optimization objective, which represents the “weighted total influence of edge prediction residuals in the non-converged subgraph”. The second term To prevent slack variables and Approaching 0, all constraints constructed from formulas (A1) to (A4) are strictly feasible. γ is the barrier factor. When the value of γ is small, the problem is closer to the original optimization problem. When the value of γ is large, the problem is closer to satisfying the constraints but farther from the original goal.

[0199] On the basis of satisfying the constraints, find a node consistency convergence state (x i ), neighbor voting connection (y ij ) and the minimum residual magnitude (z i ) so that the impact of the overall unconverged subgraph is minimized and all constraints have "internally feasible".

[0200] The constructed interior point method barrier function target model is to output the predicted state x by solving the consensus residual constraint and the consistency judgment of the subgraph structure. i , and these state variables guide whether to activate synchronization later.

[0201] In step S606, the initial barrier factor γ is set to 1. For all nodes in the converged state in the consistency graph, their boundary variables do not need communication synchronization in the next time step t+1, and are directly self-predicted by the edge cache to participate in the simulation solution; otherwise, only the unconverged connected subgraph area in the consistency graph is activated to trigger boundary synchronization and communication update. The barrier factor is set to γ ​​in each round of iteration. t+1 =η·γ t , γ t+1 and γ t are the barrier factors of the t+1th time step (i.e., the t+1th round iteration, or the t+1th generation iteration) and the tth time step, respectively. η is the decay factor, η = 0.5. The minimum iteration termination condition for boundary synchronization and communication update is that the barrier factor γ is less than 10 -4 or when the convergence tolerance is met.

[0202] Figure 16The cache content IDs of different edge nodes Edge1, Edge2, and Edge3 using the local priority strategy, the edge node caching mechanism disclosed in this disclosure that constructs a coupled boundary prediction consistency graph and performs graph-driven cache reuse scheduling (specifically steps S602-S606), and the centralized strategy are shown. The specific cache content IDs are shown in Table 1.

[0203] Table 1

[0204]

[0205] By combining Figure 16 As shown in Table 1, for the local optimal strategy, the content IDs of Edge1 nodes: [2, 19, 22, 24] and Edge3 nodes: [2, 19, 23, 25] are partially repeated, the content IDs of Edge1 nodes and Edge2 nodes have no overlap at all, and the content IDs of Edge2 nodes and Edge3 nodes also have no overlap at all.

[0206] For the centralized strategy, the content IDs of Edge1, Edge2, and Edge3 nodes completely overlap, without any collaboration. Each node only caches its own high-frequency content, resulting in high content overlap, high redundancy, and low total cache coverage.

[0207] For the disclosed method, the content IDs of Edge1 and Edge2 nodes have no overlap at all, the content IDs of Edge1 and Edge3 nodes only overlap by "19", and the content IDs of Edge2 and Edge3 nodes have no overlap at all, indicating that the content is divided and shared according to consistent subgraphs. At the same time, after the inconsistent subgraphs are discretized, the unconverged connected subgraph areas are activated to reduce invalid repetitions and retain effective overlaps. The redundancy is moderate but intelligent and selective.

[0208] The cache content ID is the content stored in the boundary cache area of ​​each edge node. It is the predicted boundary variable value, subdomain status data, or boundary interaction content during the simulation process. If multiple edge nodes have the same cache content, redundancy will occur. In particular, if the cache is repeated between nodes that do not need to be synchronized, it will be invalid redundancy. Each content ID represents a cache entity of specific physical boundary prediction data (such as temperature, potential, displacement, electric field, or flow rate parameters) at one end. It is used for cross-subdomain interaction and synchronization to reduce recalculation or transmission. If the cache content is distributed discretely (with little overlap), the overall content coverage is high and redundancy is low. If the cache content is highly repeated, it indicates a lack of coordination or sharing mechanism between adjacent nodes.

[0209] Figure 16 In the simulation, the central server CS node is responsible for maintaining all content IDs in the simulation task;

[0210] Edge1 to Edge3 represent three edge nodes, each deploying local edge cache content. The BS site represents a regional coordination base station, which serves as the edge cache center or transit coordination node for multiple edge cache nodes Edge1 to Edge3 and is used to relay synchronization requirements between edge nodes. The red solid line indicates that there is a communication requirement in an unconverged edge state between the Edge1 to Edge3 edge nodes, requiring synchronization. The blue dotted line represents the backup path for transmitting backup cache data from the BS site to the Edge1 to Edge3 edge nodes. This structure enables intelligent collaborative reuse of edge cache content and reduces communication overhead on the edge side.

[0211] For network data transmission consisting of the CS site (central server), BS site, and Edge1, Edge2, and Edge3 nodes, the redundancy rate is calculated as follows:

[0212]

[0213] For the local optimal strategy, there are 4×3=12 content IDs in total, and the total number of contents (the number of non-overlapping content IDs after deduplication) is 10. Then the local optimal strategy

[0214] For the 12 content IDs in some embodiments of the disclosed method, after deduplication of "19" where only the content IDs of Edge1 and Edge3 nodes overlap, the number of non-overlapping content IDs is 11, and the redundancy rate of the disclosed method is 8.3%;

[0215] For the centralized strategy: all nodes cache the same content, and the number of non-overlapping content IDs is 4, then the centralized strategy

[0216] For the effective content coverage indicator, the calculation formula is:

[0217] If the content library of the CS site has 30 items, then the coverage rate for the local optimal strategy is: 10 / 30=33.3%, the coverage rate for the disclosed method is: 11 / 30=36.7%, and the coverage rate for the centralized strategy is: 4 / 30=13.3%.

[0218] The redundancy rate and coverage of the three methods are summarized in Table 2.

[0219] Table 2

[0220] Strategy Redundancy rate Effective content coverage Local optimization strategy 16.7％ 33.3％ Centralized strategy 66.7％ 13.3％ The disclosed method 8.3％ 36.7％

[0221] The edge node caching mechanism proposed in this disclosure, which builds a coupled boundary prediction consistency graph and performs graph-driven cache reuse scheduling, demonstrates the following advantages, compared to existing local optimal and centralized strategies, in the same simulated network structure: effectively reducing cache redundancy and avoiding duplicate caching. Using this embodiment of the disclosure, each edge node only performs boundary cache communication and synchronization in areas where predictions are inconsistent, and reuses local cache in areas where predictions converge.

[0222] Figure 16 The results in Tables 1 and 2 show that in the traditional centralized strategy, multiple nodes repeatedly cache exactly the same content (for example, nodes Edge1 to Edge3 all cache [3, 19, 22, 29]), resulting in a redundancy rate of up to 66.7% in the edge cache resources. However, the disclosed method drives node collaboration through the graph structure, sharing the cache only in inconsistent subgraphs, reducing the redundancy rate to 8.3%, and covering more unique content (with the highest coverage of effective content, reaching 36.7%).

[0223] Furthermore, the present disclosure can also improve content coverage and enhance the flexibility of reuse scheduling. Unlike centralized strategies that only cache fixed centrally distributed content, the present disclosure utilizes a dynamic subgraph partitioning mechanism based on the predicted consistency graph, enabling each edge node to participate in caching tasks based on its own state and neighborhood residuals. Therefore, the method of the present disclosure can improve the uniqueness and distribution diversity of boundary cache content, achieving maximum boundary state coverage within the same cache space constraints. Furthermore, when multiple simulation subdomains are collaboratively computed, the method effectively supports cache reuse across more sub-physical domain boundaries, improving simulation coupling efficiency.

[0224] The present disclosure implements a local consistency consensus voting process by introducing an interior point method barrier function target model constrained by formulas (A1) to (A4) in step S604, and obtains a global decision on whether the node is in a predicted convergence state, ensuring that communication updates are activated only in areas with inconsistent predictions, and that boundary nodes in consistent subgraphs all use cached prediction values ​​instead of communication; compared to centralized solutions that require frequent full-graph synchronization, the method of this embodiment of the present disclosure can avoid approximately 60% or more of unnecessary communication path activations, especially in high-bandwidth load or cluster environments, significantly reducing delays and computing overheads. Therefore, the present disclosure also achieves controllable communication costs and significantly reduces the number of global synchronizations.

[0225] Figure 17 A flowchart of a method 700 for constructing an adaptive scheduling diagram and dynamically optimizing simulation resource allocation and solution strategies based on real-time status monitoring results of simulation tasks according to an exemplary implementation of the present disclosure is schematically shown.

[0226] In step S702, for each node v i Defining Embedding Features Form an adaptive scheduling graph;

[0227]

[0228] Integrates the i-th node v i As the state variable μ of the i-th physical domain where the vertex is located phys,i (such as temperature, displacement, electric potential, etc.), current coupling centrality C i And the communication complexity Comm i , network size, boundary variable history residuals, residual diagrams, regional boundary tension / disturbance intensity (such as heat flow jumps, current density mutations, etc.), state history change rate, local transmission rate or unit boundary communication bandwidth weight factor, etc. The "..." in the embedded feature represents the state variables that can also be included. μ phys,i Entered as part of the graph embedding vector It is the initial input feature for graph neural network (GNN) scheduling optimization, information propagation analysis and bottleneck identification.

[0229] In step S704, based on the GNN information propagation model, multiple rounds of information aggregation between different physical domains on the adaptive scheduling graph are simulated:

[0230]

[0231] For two different nodes v i and v j , both of which satisfy v i ∈V,v j ∈V, the nodes adjacent to the i-th node in the constructed adaptive scheduling graph, and all nodes v adjacent to the i-th node j The set of N is defined as F (v j ): j-th node v j The graph distance is defined as d F (v i ,v j ), F is the scheduling graph of information propagation between the i-th node and all its adjacent nodes; Represents the i-th node v in the l-th layer of the graph neural network i The aggregation result is obtained by the function f l (·) For the i-th node v in the l-1th layer i itself and the l-1th layer aggregation result of its adjacent node set To express the obtained.

[0232] In a graph neural network with L layers, node v i The final embedding It can obtain structural information within a range of up to L hops (“hop” is the distance (i.e., the number of edges) from a certain node to the node that can be reached through layer-by-layer propagation of edges). i Starting from the expansion aggregation step, we can get a computational graph with as the root, which is a tree of depth L, representing an L-hop neighborhood, where the children of any node in the tree are the adjacent nodes.

[0233] Figure 18 A flow chart of a method 800 for identifying information transmission bottleneck areas based on an information contraction criterion according to an exemplary implementation of the present disclosure is schematically shown. In some embodiments of the present disclosure, method 800 may be a method flow specifically included in step S706. In some embodiments, step S706 specifically includes steps S802-S806. In step 806, the information contraction (information bottleneck) criterion is used to calculate the information loss rate of the key node / area:

[0234] The i-th node v i After multiple rounds of iterations, the amount of global coupling information obtained should be gradually shrunk by introducing the information shrinkage coefficient η L (E) is used to measure the attenuation of information arrival in the neighborhood (i.e., the physical domain with multiple adjacent nodes as vertices). E is the information propagation path in the Markov chain adaptive scheduling graph from U→A→B, where U is the global initial variable or key coupling physical quantity of the simulation system. B is the i-th node v i The ith physical domain where A is located is the ith node v i One of the multiple physical domains corresponding to the adjacent multiple nodes (with the jth node v j For example, U→A→B represents the physical information flow path of the simulation information from the global state through region A to region B. Therefore, channel E represents all information propagation paths from region A to region B on the scheduling graph F corresponding to the scheduling graph. The bottleneck region S is screened by changing the entropy or mutual information embedded in each node at different levels.

[0235] Define the mutual information between the global variable U and the middle area A of the channel of the information propagation network as I(U;A), and the mutual information between the global variable U and the target area B (for example, the target node v i The mutual information of is I(U; B), if the ratio of I(U; B) to I(U; A) is less than the information shrinkage coefficient η L (E), it indicates that the path or regional structure U→A→B may have channel interference, over-compression, too many layers, and other problems; at this time, the bottleneck area S of the information propagation channel E is identified, and its associated channel E is recorded as the path to be reconstructed.

[0236] That is, the condition is when When , it is determined that there is a bottleneck area S in the information propagation channel E.

[0237]

[0238] Figure 19 (a) Schematic diagram of the multi-hop neighborhood aggregation structure based on graph embedding features during information propagation, where each node in the figure represents a computing unit in the simulation task, the arrow indicates the propagation direction of the state information, and the nested subgraph structure represents the dense coupling characteristics of the internal structure of the bottleneck area; The aggregated state information of node B is output to the top layer in sequence, reflecting the path process of the graph neural network to aggregate the global state layer by layer under the advancement of embedding depth. L In the calculation formula of (E), sup[] represents the supremum function. are the nth possible state distribution and the mth possible state distribution of area A respectively, (i.e. ) is the path distribution from area A to area B under the nth possible state distribution, (i.e. ) is the path distribution of region A mapping to region B under the mth possible state distribution, n = 1, 2, ..., M; m = 1, 2, ..., M. η L (E)∈[0,1],η L The smaller (E), the more serious the information loss is during the information transmission from area A to area B, i.e., bottleneck or over-compression. is the information separation measure of the two possible states of region A, That is to say It is a measure of information separation under the above two path distributions, and its calculation formula is the KL divergence calculation rule and formula.

[0239] for and The input state space has N state variables, so the probability of the nth and mth possible state distributions in area A is as follows:

[0240] P(a ni )for The state variables of the i-th physical domain in the input state space;

[0241] P(a mi )for Enter the state variables of the ith physical domain in the state space. That is, a ni with ami Represents different distribution values ​​under the same physical state variable (i.e., the variables correspond one to one).

[0242] The result of mapping different possible state distributions of area A to area B is the set of N state variables in area B:

[0243] The probability matrix of the mapping of channel E in the Markov chain adaptive scheduling graph of U→A→B is:

[0244]

[0245] The distribution probability of two different possible states of area A being transmitted to area B through channel E Think of it as a matrix-vector multiplication of an N-dimensional vector. Define the two output edge distributions that map different possible state distributions in area A to area B.

[0246]

[0247] Information separation metrics under two path distributions

[0248]

[0249] The following is a specific calculation method using N=2: but

[0250]

[0251] for The calculation rules are the same as above. Then the probability matrix of the mapping of channel E in the Markov chain adaptive scheduling graph of U→A→B is:

[0252]

[0253] For example, when So Then the numerator calculation result of the information shrinkage coefficient is:

[0254]

[0255] The denominator of the information shrinkage coefficient is calculated as:

[0256] and then

[0257] Mutual information can be calculated using probability distribution (embedding entropy) or entropy-based methods in graph neural networks, or with the help of deep embedding vector mutual information estimators (such as MINE and JS-MI). Taking the entropy-based method in graph neural networks as an example:

[0258] I(U;A)=H(A)-H(A|U);

[0259] H(A) is the total uncertainty (the degree of information confusion) of region A in the graph neural network:

[0260] H(A)=-∑ a∈A P(a)logP(a); P(a) is the probability of each state vector a appearing in region A (which can be obtained by statistically analyzing historical samples). H(A|U) represents the remaining uncertainty entropy value of region A after the global variable U is known.

[0261]

[0262] Figure 19 (b) shows that for the bottleneck area S obtained by screening, the internal node connectivity is poor, resulting in information transmission attenuation; therefore, it is necessary to perform structural fitting analysis and graph reconstruction on the bottleneck area.

[0263] Figure 20 A flow chart schematically illustrates a method 900 for performing structural fitting analysis and atlas reconstruction on a bottleneck region according to an exemplary implementation of the present disclosure. In some embodiments of the present disclosure, the method 900 may be used as a method flow specifically included in step S708, i.e., in some embodiments, step S708 specifically includes steps S902-S906.

[0264] In step S902, a regular graph g is selected or constructed for the bottleneck region S. S :g S =(V S ,E S ) is a p-regular graph, that is, a regular graph g S Each node in is connected to p edges. For a regular graph, its adjacency matrix is ​​Ad S , Ad S The element Ad in row i and column j si,sj as follows:

[0265]

[0266] Ad si,sj =1, indicating that the regular graph g S The i-th node v in si Connected to the jth node; Ad si,sj =0, indicating that the regular graph g SThe i-th node v in si Not connected to the jth node;

[0267] Define the adjacency matrix Ad S The eigenvalue spectrum of Ad S e=λ n ·e;

[0268] Among them, e is the regular graph g S The node connection feature vector in ,λ n is a regular graph g S The nth eigenvalue in ;

[0269] Assume that the order of eigenvalues ​​is λ1≥λ2≥…≥λ n ; For the regular graph g S For example, the main eigenvalue of its adjacency matrix is ​​λ1; λ2 is the second largest eigenvalue of the adjacency matrix.

[0270] Then in step S904, the regular graph g S The local spectral gap Δ spec (g S )=d s -λ2,d s is a regular graph g S the average degree of the nodes in (usually computed using linear algebra tools such as numpy.linalg.eigvals(Ad_S)),

[0271] deg(v si ) represents the computing node v i The number of edges connected in the graph, that is, the number of adjacent nodes that determine the node. |S| is the total number of nodes in the regular graph.

[0272] If deg(v si )=d s (i.e. each node v si The number of connected edges in the graph is equal), then the original graph S is approximately regular, and its adjacency matrix Ad can be directly extracted S .

[0273] If S is not regular, then the regular graph g S Replaced by a structurally fitting d-regular graph prototype (such as Figure 21 ring-of cliques, equidetermined random graphs, etc.)

[0274] λ2 reflects the tightness of node connections inside the bottleneck region S. If the regular graph g S The local spectral gap Δ spec (g S) is significantly smaller than the average spectrum gap of the entire image, it indicates that there is a structural bottleneck in this area, and it is necessary to trigger the random local edge flipping and edge optimization reconstruction.

[0275] In step S906, for the bottleneck region S identified in step S706, if its local spectrum gap Δ spec (g S ) is significantly smaller than the average spectrum gap Δ spec (G) indicates that the internal connections in the region are sparse (i.e., the nodes are not tightly connected) and the structure is severely contracted. Therefore, it is necessary to reconstruct the edges in this step. The specific steps for reconstructing the connections are as follows:

[0276] First, the local spectral gap Δ obtained by screening spec (g S ) is significantly smaller than the average spectral gap of the full graph. S Select the node pairs (u,v) with a small number of connected edges (degree) or a weak change in the embedded information feature as potential reconstruction targets. S In some embodiments, the protein is obtained by screening in step S706.

[0277] Then, for the candidate edge (u, v), select adjacent nodes i∈N(u) and j∈N(v) from their respective neighbor node sets, satisfying one or more of the following conditions:

[0278] 1) Satisfy the symmetry of the coupling relationship (such as matching of temperature and potential boundary conditions);

[0279] 2) Satisfy topological constraints such as edge length and network hop count (e.g., within the maximum hop count L). If so, mark the edge pair (i,u), (j,v) as a candidate for reconstruction. 3. Edge flip operation: Delete the edges (i,u), (j,v) from the edge set E in the original graph and replace them with the new edges (i,v), (j,u). The merged results form a new edge set, and the new graph structure is constructed:

[0280] G′=(E\{(i,u),(j,v)})∪{(i,v),(j,u)};

[0281] Among them, {(i,u),(j,v)} are the old connection edges to be removed; {(i,v),(j,u) are the newly generated reconstructed edges; ∪ is the union; G′ represents the optimized subgraph generated after the edge flipping operation.

[0282] Through the above edge reconstruction connection process, ensure that the bottleneck area g S The information flow path in the network is optimized, thereby improving the connection density and spectrum gap index of the region, and improving the graph propagation efficiency and resource coordination ability in the entire simulation task.

[0283] about Figure 21 ring-of cilques, such as Figure 21 (b) for Figure 21 The original graph (non-regular) of (a) is used to identify subgroup structures or neighbor communities, and the nodes in the original graph are divided into communities or blocks (for example, the non-regular nodes in the original graph are constructed into two subgroups: {0, 1, 2, 3} and {4, 5, 6, 7}). The grouping criteria can be based on degree, adjacency density or physical proximity. Figure 21 The grouping criterion for (b) is that the final regularity is required to be 4. Then construct a clique (complete graph) for each subgroup, and for each group of nodes, construct a complete graph (each node is connected to all other nodes in the group) to form a clique. Assuming that each group has k nodes, the clique is a k-clique. Figure 21 Each subgroup in (b) has 4 nodes. Connect the cliques in a ring to form a "ring", select a "representative node" from each clique and connect it to the representative node of the next clique to form a ring structure, corresponding to Figure 21 In example (b), two 4-cliques are connected in sequence through bridges of their respective subgroup nodes (i.e., node 0 is connected to node 4, node 1 is connected to node 5, node 2 is connected to node 6, and node 3 is connected to node 7), forming a ring-of-cliques, thereby ensuring that the degree of all nodes in the graph remains consistent (e.g., d = 3), thereby achieving d-regularity.

[0284] Figure 21 (c) shows the Figure 21 The original graph of (a) is constructed (i.e. fitted) by constructing an equal-degree random graph. Figure 21 The original graph of (a) retains the original graph node set and replaces all nodes in the original graph (such as Figure 21 (a) The nodes numbered 0-7 are retained and not clustered or grouped. Specify the target regularity d - Select the target regularity d, such as d = 3, which means that each node is connected to 3 edges. Use the random regular graph generation algorithm (such as configuration model or networkx.random_regular_graph(d,n)) to generate an equal-degree random graph of the relevant nodes of the original graph, ensuring that all nodes have degree d, and there are no duplicate edges and no self-loops, thus forming a Figure 21 (c) is an equal-degree random graph of red nodes, each node has 3 adjacent edges, Figure 21 (c) The overall graph is a 3-regular graph, the edge connections are completely random, and the structure of the original graph is not preserved.

[0285] Figure 21 The adjacency matrix of the original graph in (a) is like Figure 22 As shown, the number of rows in the matrix from left to right and each column from top to bottom are arranged in order from node 0 to node 7. If there is a connecting bridge between the two, the value of the element located at the matrix position of the element is 1, otherwise it is 0. The eigenvalues ​​of each node are shown in Table 3:

[0286] Table 3. Sorting of eigenvalues ​​of original graph

[0287] Eigenvalue number Eigenvalue <![CDATA[λ 1,a ]]> 2.3028 <![CDATA[λ 2,a ]]> 0.618 <![CDATA[λ 3,a ]]> 0.0 <![CDATA[λ 4,a ]]> -1.3028 <![CDATA[λ 5,a ]]> -1.618

[0288] The average degree of nodes in the original graph is:

[0289]

[0290] Then the global spectrum gap of the bottleneck region Sa based on the original graph is Δ spec (G a )=d S,a -λ 2,a =2.0-0.618=1.328. This value is subsequently used to determine whether the spectral gap of a local graph (such as a ring-of-cliques regular graph or an iso-degree random graph) is "significantly smaller" than the average spectral gap of the entire graph, thereby identifying whether the bottleneck region still exists and triggering structural reconstruction.

[0291] Figure 21 The adjacency matrix of the ring-of clique regular graph in (b) is The characteristic values ​​of each node are shown in Table 4:

[0292] Table 4. Ranking of eigenvalues ​​of ring-of-Cliques regular graphs

[0293] Eigenvalue number Eigenvalue <![CDATA[λ 1,b ]]> 4.0 <![CDATA[λ 2,b ]]> 2.0 <![CDATA[λ 3,b ]]> 0.0 <![CDATA[λ 4,b ]]> -0.0 <![CDATA[λ 5,b ]]> -0.0 <![CDATA[λ 6,b ]]> -2.0 <![CDATA[λ 7,b ]]> -2.0 <![CDATA[λ 8,b ]]> -2.0

[0294] According to the above node average degree calculation formula, the graph is Figure 21 (b) The average node degree d of the ring-of-Cliques regular graph S,b is 3, the local spectrum gap Δ spec (g S,b )=d S,b -λ 2,b =4.0-2.0=2.0.

[0295] Δ spec (g S,b )>1.328, so if we follow Figure 21 The regular graph constructed in the manner of (b) does not need to be reconstructed.

[0296] Figure 21 The adjacency matrix in the equal-degree random graph of (c) is The characteristic values ​​of each node are shown in Table 5:

[0297] Table 5 Ranking of eigenvalues ​​of equal-degree random graphs

[0298] Eigenvalue number Eigenvalue <![CDATA[λ 1,b ]]> 3.0 <![CDATA[λ 2,b ]]> 1.0 <![CDATA[λ 3,b ]]> 1.0 <![CDATA[λ 4,b ]]> 1.0 <![CDATA[λ 5,b ]]> -1.0 <![CDATA[λ 6,b ]]> -1.0 <![CDATA[λ 7,b ]]> -1.0 <![CDATA[λ 8,b ]]> -3.0

[0299] According to the above node average degree calculation formula, the node average degree d of the ring-of-Cliques regular graph in (b) is S,b is 3, the local spectrum gap Δ spec (g S,b )=d S,b -λ 2,b =3.0-1.0=2.0.

[0300] Δ spec (g S,c )>1.328, so if we follow Figure 21 The regular graph constructed in the manner of (b) also does not need to be reconstructed.

[0301] The disclosed method combines the propagation principles of graph neural networks, proposes a bottleneck identification method based on information shrinkage ratio and local spectral gap, and innovatively introduces a graph structure reconstruction mechanism (such as edge flipping and regularized fitting), realizing the self-evolution of communication topology and automatic optimization of weakly connected areas, thereby improving the information flow capacity of multi-physics field systems.

[0302] The method 900 disclosed in the present invention can accurately locate and eliminate communication and information transmission bottlenecks in the system, thereby enhancing the critical path connectivity between nodes during the simulation process, greatly improving the overall data flow efficiency and robustness of the system, and providing a foundation for subsequent complex large-scale parallel simulation and intelligent scheduling optimization.

[0303] although Figure 7 、 Figure 8 、 Figure 9 、 Figure 14 、 Figure 15 、 Figure 17 、 Figure 18 and / or Figure 20 The various steps shown in the embodiment perform the corresponding functions of the performance recommendation method in a specific order, but it can be understood that Figure 7 、 Figure 8 、 Figure 9 、 Figure 14 、 Figure 15 、 Figure 17 、 Figure 18 and / or Figure 20The steps in the method 200 may be performed in a different order relative to the order shown. For example, two or more steps shown in succession may be performed simultaneously, or portions thereof may be performed simultaneously. Furthermore, certain steps in the method 200 may be omitted, as may be the case with methods 300, 400, 500, 600, 700, 800, and 900.

[0304] Additionally, other steps can be added to Figure 7 、 Figure 8 、 Figure 9 、 Figure 14 、 Figure 15 、 Figure 17 、 Figure 18 and / or Figure 20 It will be understood that all such variations are within the scope of the present disclosure.

[0305] The present invention also proposes an electronic device 10 suitable for accelerating multi-physics field coupling simulation for industrial digital twins, such as Figure 23 As shown. The electronic device includes at least one processing unit 101, at least one storage unit 102, and a computing array integrated with multiple simulation computing nodes. The processing unit 101 can be a general-purpose CPU, GPU, FPGA or other high-performance computing unit, which is used to execute core algorithms such as simulation master control process, parallel computing scheduling, state monitoring, data fusion and graph structure optimization. The storage unit 102 may include a cache, main memory (RAM) and / or non-volatile memory (such as SSD, hard disk) for storing information such as simulation models, task queues, intermediate calculation data and boundary variable cache.

[0306] The computing array integrates multiple computing nodes (such as nodes A, B, C, D, E, and F). Each node contains multiple computing subunits (such as T1–T8, U1–U8, V1–V8, W1–W8, H1–H8, and R1–R8), which are used to execute local solution tasks in their respective subdomains in parallel. Multi-threaded / multi-core parallel computing is supported within the node, enabling refined subdomain decomposition and collaborative solution. Each node can correspond to a physical domain or simulation subdomain, enabling task decomposition and parallel scheduling between subdomains.

[0307] In some embodiments, the electronic device provided by the present disclosure includes Figure 10-12 Multiple computing nodes, Figure 10-12 The computing sub-units between multiple computing nodes in the Figure 13In the computing array shown, simulation nodes are interconnected via high-speed buses, dedicated communication interfaces, or a network-on-chip (NoC). These interfaces can utilize PCIe, Ethernet, InfiniBand, or high-speed interconnect protocols within the SoC to ensure low-latency, high-bandwidth transmission of state variables and boundary data between nodes. Furthermore, the electronic devices can reserve external expansion interfaces, such as USB, Ethernet, and serial ports, for data exchange and real-time control with external sensors, hosts, or industrial control systems.

[0308] The electronic device 10 disclosed in the present invention can adaptively support different node communication connection topologies according to simulation requirements and task distribution. Specifically, each computing node in the computing array can not only solve local tasks of a single node, but also flexibly construct various types of inter-node data communication paths including central cross, regular grid, redundant multi-ring, etc. through high-speed interconnection mechanism, thereby achieving Figure 10 – Figure 12 The diverse task scheduling and regional coupling relationships shown in Figure 1 are shown. The electronic device 10 can automatically configure inter-node communication connections based on the structural characteristics, coupling strength, or scheduling strategy of the multi-physics coupling model, dynamically optimizing overall simulation efficiency and load balancing. This architecture provides a unified, scalable hardware and system foundation for large-scale parallel simulation and efficient processing of heterogeneous coupling problems.

[0309] Figure 23 The storage unit 102 is pre-installed with computer-readable program instructions. When loaded and executed by the processing unit 101, it drives the above-mentioned nodes to implement the multi-physics field coupling simulation acceleration steps of the above-mentioned embodiment, including: sub-domain division, graph construction, scheduling optimization, cache reuse, adaptive bottleneck identification and structure reconstruction and other functions.

[0310] In summary, this electronic device provides hardware support and expansion capabilities for the industrial digital twin system with its highly integrated multi-node simulation array, flexible interface architecture and efficient data transmission capabilities, ensuring the efficient and reliable operation of multi-physics field coupling simulation tasks.

[0311] Those skilled in the art will appreciate that the method steps described herein are not limited to the order exemplarily shown in the drawings, but may be executed in any other feasible order.

[0312] The above description is provided to enable any person skilled in the art to make or use the present disclosure. Various modifications of the present disclosure will be readily apparent to those skilled in the art, and the general principles defined herein may be applied to other variations without departing from the spirit and scope of the present disclosure. Therefore, the present disclosure is not limited to the examples and designs described herein, but is intended to be consistent with the widest scope of the principles and novel features disclosed herein.

Claims

1. A multi-physics field coupling simulation method for industrial digital twins, characterized by: The following steps are involved: Constructing a coupled simulation model based on the actual physical field data of industrial equipment or systems, dividing the different physical domains in the coupled model into sub-models, and establishing the solution area and boundary conditions for each physical domain; The coupled model is divided into regions using a domain decomposition method to construct a corresponding region-level parallel computing graph, where each region serves as a computing node and the shared boundaries between nodes serve as connecting edges, forming a topological structure for synchronous solution of multiple physical fields; Based on the parallel graph, multi-level scheduling and computing task mapping are performed; During the simulation process, an edge node cache mechanism is used to store key state variables of each sub-model, and the cached variables are reused in the cross-domain solution synchronization phase to reduce repeated calculations and communication delays; Based on the real-time status monitoring results of simulation tasks, an adaptive scheduling diagram is constructed to dynamically optimize simulation resource allocation and solution strategies to achieve accelerated operation of multi-physics field modeling and simulation in industrial digital twin systems.

2. The method according to claim 1, characterized in that The coupled model is divided into regions to form a topological structure for simultaneous solution of multiple physical fields, including: Construct a variable dependency graph, where vertices represent state variables or physical subdomains, edges represent coupling driving relationships between variables, and edge weights are calculated based on the partial derivatives of the residuals of the control equations with respect to the coupled variables. The coupling centrality of each variable is calculated based on the edge weights in the graph, reflecting the strength of its association with other variables. The priority of computing nodes is determined based on the coupling centrality results, and highly coupled nodes are allocated to areas with stronger computing power or higher communication bandwidth. Establish a scheduling objective function to minimize the cross-region communication cost between highly coupled nodes and optimize the parallel computing topology.

3. The method according to claim 1, characterized in that The multi-level scheduling and computing task mapping includes: At the first level, a load balancing strategy is used to build a computation priority scheduling queue based on the computational complexity and communication dependency of the physical fields in each region; At the second level, a graph partitioning algorithm is used to map regional computing tasks to processing units or GPU clusters on a multi-core computing platform to form a computing resource assignment table. At the third level, based on the time stepping and synchronous triggering mechanism, the iterative process of each area is integrated with tasks, compressed with communications, and scheduled with synchronous solutions to achieve collaborative evolutionary simulation of multiple physical domains.

4. The method according to claim 3, characterized in that The first level includes: Obtain the coupling centrality and communication complexity indicators of each regional node; Calculate the comprehensive scheduling weight of each node to balance computing load and communication intensity; Sort all regional nodes according to their weight values ​​and build a scheduling priority queue; Each computing task node in the scheduling queue is submitted to the task scheduler in sequence and waits to be allocated to the computing resource unit according to priority.

5. The method according to claim 3, characterized in that The third level includes: Execute the solution of local coupling control equations in each physical area and complete single-step iterative calculations; compress and transmit boundary variables between areas to reduce the amount of communication data; Determine whether the inter-regional coupling residuals meet the synchronous advancement conditions and control the time step consistency; If convergence has not occurred, a new round of iterative scheduling is triggered until the co-evolution requirements are met.

6. The method according to claim 1, characterized in that The edge node caching mechanism also includes building a coupled boundary prediction consistency graph and performing graph-driven cache reuse scheduling, specifically including: In each simulation iteration cycle, a cross-region prediction consistency graph structure is constructed based on the coupling boundary relationship between sub-models. Each node in the graph corresponds to a regional boundary cache, and the edge weight represents the correlation of the deviation change trend between the prediction values ​​of two boundary caches. For each node, based on its historical cache evolution sequence and the synchronization residual function between adjacent nodes, it iteratively updates its consistency state flag and performs a local consistency consensus voting process to obtain a global decision on whether the node is in the predicted convergence state; For all nodes in the consistency graph that are in a converged state, their boundary variables do not require communication synchronization in the next time step, and are directly self-predicted by the edge cache to participate in the simulation solution; otherwise, only the unconverged connected subgraph area in the consistency graph is activated to trigger boundary synchronization and communication update.

7. The method according to claim 1, characterized in that The steps of constructing and optimizing the adaptive scheduling graph include: Construct a graph structure that integrates physical states and coupling indicators for scheduling optimization; Perform multiple rounds of information propagation through graph neural networks to extract key path structural features; Identify information transfer bottleneck areas based on information contraction criteria; Perform structural fitting analysis and graph reconstruction on bottleneck areas to optimize node connectivity and resource scheduling efficiency.

8. The method according to claim 1, characterized in that The identifying of the information transmission bottleneck area based on the information contraction criterion includes: Based on the graph neural network model, multi-level neighborhood information aggregation operations are performed on the adaptive scheduling graph; Gradually expand the node structure perception range to obtain the state expression within the multi-hop neighborhood; Based on the information contraction theory, the information loss rate of key nodes is calculated to determine whether there is a bottleneck area in the information propagation channel.

9. The method according to claim 6, characterized in that The performing of structural fitting analysis and atlas reconstruction on the bottleneck region includes: Select or fit a regular graph structure for the bottleneck area obtained by screening, extract its adjacency matrix and spectral features; calculate the local spectral gap index of the regular graph to measure the tightness of node connections and structural balance; When the local spectral gap index is significantly lower than the average level of the entire graph, it is judged that the nodes in this area are not tightly connected and subsequent edge structure optimization processing is required.

10. An electronic device comprising at least one processing unit and at least one storage unit, wherein: The storage unit stores a computer program, and when the computer program is executed by the processing unit, the processing unit is caused to perform the steps of the method according to any one of claims 1 to 9.

Images