A method for improving the accuracy of multiphysics coupled simulation solutions

Abstract

This invention relates to the fields of computer-aided engineering and high-performance computing, and discloses a method for improving the accuracy of multiphysics coupled simulation solutions. This method first divides the computational domain based on the underlying hardware architecture of the computing device to optimize data locality. During the solution process, the accumulated energy residual at the coupling interface is monitored in real time. By comparing its rate of change with the pre-trained instability critical slope, a proactive warning of computational divergence is achieved. Once a warning is issued, the system transforms the residual into a physically meaningful additional nodal force vector, which is then smoothly applied back to the target field solver in the form of a dynamic source term, combined with spatial weights, to achieve physically consistent error correction. This invention improves computational efficiency through hardware-aware partitioning, saves computational resources through a proactive warning mechanism, and enhances the accuracy and reliability of simulation results through a physically consistent energy correction scheme.

Description

Technical Field

[0001] This invention relates to the fields of computer-aided engineering and high-performance computing, specifically to a method for improving the accuracy of multiphysics coupled simulation solutions. Background Technology

[0002] Multiphysics coupled simulation is a key tool in modern engineering design and scientific research. It simulates the interactions in complex systems by having solvers from different physical domains perform collaborative calculations. In traditional coupled solution strategies, each physical solver proceeds within an independent computational framework and exchanges data through a pre-defined coupling interface. For example, the fluid solver transmits pressure to the structural solver, which in turn returns displacement information.

[0003] However, existing technologies face inherent challenges in dealing with this coupling problem. On the one hand, when data is transferred at the coupling interface, it is difficult to ensure strict energy conservation due to grid mismatch between different physical fields, inconsistent time steps, and inherent errors in interpolation algorithms. This small energy imbalance will accumulate continuously over thousands of time steps, eventually leading to the generation or dissipation of non-physical energy. This not only reduces the accuracy of simulation results, but in severe cases, it can also cause numerical instability, causing the entire calculation process to be interrupted prematurely.

[0004] On the other hand, as the scale of computing increases, these simulation tasks increasingly rely on massively parallel computing devices. When performing domain decomposition, existing simulation software only considers the load balancing of the physical model, while ignoring the underlying hardware topology of the computing devices, such as non-uniform memory access architecture. This mismatch between hardware and software leads to a large number of cross-node remote memory accesses, and the resulting communication latency becomes a bottleneck restricting the improvement of overall computing efficiency.

[0005] Furthermore, when simulations exhibit instability, existing monitoring and correction mechanisms are rather passive. They typically only trigger intervention when the residual of a physical quantity exceeds a fixed, empirically set threshold. By this time, the calculation may have deviated significantly from the correct solution. The corrective measures taken, such as simply applying a reverse force or resetting the time step, are rather rigid and introduce new numerical oscillations into the model, interfering with the smooth convergence process of the solution. Therefore, designing a method that can synergistically improve computational efficiency, proactively maintain solution stability, and perform error correction in a physically consistent manner is a technical problem that urgently needs to be solved in this field. Summary of the Invention

[0006] To address the shortcomings of existing technologies, this invention provides a method for improving the accuracy of multiphysics coupled simulation solutions. This method solves the problem that existing technologies in multiphysics coupled simulations are unable to simultaneously improve simulation accuracy, computational efficiency, and numerical stability due to energy non-conservation, mismatch between computational tasks and underlying hardware, and lack of physical consistency in error correction.

[0007] To achieve the above objectives, the present invention provides the following technical solution: a method for improving the accuracy of multiphysics coupled simulation solutions, wherein the method is configured to be executed in a computing device, the computing device including a processor, a memory, and a communication bus, the processor establishing a low-level data transmission channel with the memory through the communication bus, and the method comprising the following steps: S1, the global coupling interface is divided into a local energy conservation domain according to the hardware topology boundary of the computing device, and a spatial distance attenuation weight matrix is ​​generated. S2, extract the numerical stiffness feature norm and the spatial gradient extrema of the feature physical quantity of the local energy conservation domain to calculate and generate the tolerance threshold. S3 controls the source field solver and the target field solver to advance asynchronously and performs spatial interpolation mapping at the coupling communication point; S4, calculate the difference between the virtual work output by the source field and the virtual work received by the target field and accumulate it as the cumulative energy residual; S5, determine the instability risk based on the accumulated energy residual and trigger an interrupt signal; S6, combining the spatial distance attenuation weight matrix, the accumulated energy residual is transformed into an additional nodal force vector, which is then superimposed on the external source term of the target field solver.

[0008] By adopting the above technical solution, this method combines the underlying hardware topology of the computing device, local physical field characteristics, and global energy conservation constraints. It optimizes data locality through hardware topology-aware region partitioning, providing a foundation for energy calculation and correction. The dynamically generated tolerance threshold is correlated with local physical characteristics, making error monitoring adaptive. This method constructs a closed-loop energy correction process. By monitoring the energy residual at the coupling interface, it intervenes before determining the risk of instability. The accumulated energy error is fed back to the target field solver in the form of additional nodal forces, thereby improving the solution accuracy of multiphysics coupling simulation while ensuring computational stability.

[0009] Preferably, in S1, the step of dividing the global coupling interface into a local energy-conserving domain according to the underlying hardware topology boundary of the computing device includes: Extract non-uniform memory access architecture configuration information, parse and extract the affinity distance matrix between the processor's internal processing core and different physical addressing segments of the memory, and determine the physical node boundary at the bottom layer of the computing device based on the affinity distance matrix; The multi-level graph partitioning algorithm is called to divide the global coupling interface into multiple non-overlapping and spliced ​​local energy conservation domains, and the grid node coordinates and material property data contained in each local energy conservation domain are allocated and resided in the corresponding local memory physical address segment.

[0010] By adopting the above technical solutions, cross-domain memory access latency is reduced and data exchange efficiency is improved. At the same time, the computational load balancing of local energy-conserving domains is achieved through graph partitioning algorithms, providing support for the subsequent calculation of energy residuals in each domain.

[0011] Further, in S1, the step of generating the spatial distance attenuation weight matrix includes: Collect historical perturbation propagation datasets of the global coupling interface under different grid sizes, extract the pre-obtained optimal reference weight distribution as a supervision label, construct a loss function to iteratively update the spatial decay feature length constant, and train a mathematical model of boundary distance decay kernel function. Traverse any discrete grid node within the local energy conservation domain and calculate the shortest Euclidean distance scalar to the cross-domain boundary of the corresponding local energy conservation domain. Substitute this scalar into the boundary distance decay kernel function mathematical model to generate the corresponding dimensionless weight value. Construct the spatial distance decay weight matrix using all dimensionless weight values ​​as diagonal elements.

[0012] By adopting the above technical solution, the solution uses historical data and machine learning methods to generate a boundary distance decay kernel function, so that the weight distribution reflects the propagation law of physical disturbance, thereby distributing the energy residual to the nodes near the interface according to the law, and avoiding abrupt changes in the correction force.

[0013] Preferably, in S2, the step of extracting the numerical stiffness feature norm of the local energy-conserving domain and the spatial gradient extrema of the feature physical quantities to calculate the tolerance threshold includes: Extract the system Jacobian matrix or stiffness matrix generated after discretizing the local energy conservation domain control equations, and perform Frobenius norm calculation to generate the numerical stiffness characteristic norm. In the multiphysics simulation variables, a preset master characteristic physical quantity is locked, the gradient vector of the physical quantity of the discrete grid node in three-dimensional space is calculated, and the gradient norm with the largest value inside the local energy conservation domain is selected as the spatial gradient extremum of the characteristic physical quantity. The numerical stiffness feature norm is multiplied by the spatial gradient extremum of the feature physical quantity and substituted into the empirical scaling constant obtained through data-driven training to unify the dimensions, thereby generating the tolerance threshold for local numerical stability.

[0014] By adopting the above technical solution, a dynamic tolerance threshold related to the local physical state is generated. This threshold has a larger value in regions with gentle physical fields and a smaller value in regions with drastic gradient changes, thereby achieving adaptive adjustment of monitoring sensitivity.

[0015] Furthermore, in S3, the step of asynchronously advancing the source field solver and the target field solver and performing spatial interpolation mapping at the coupling communication point includes: Extract the minimum grid size and medium wave velocity of the source field solver and the target field solver respectively, and configure the source field time step and the target field time step to be an integer fraction of the preset macroscopic coupling communication period according to the preset Courant-Friedrich-Levy stability condition. Monitor the local advance clocks of the source field solver and the target field solver, and when both are extrapolated to the corresponding coupling communication point, trigger a synchronization barrier command to pause the time extrapolation process and enter the space data exchange stage. A spatial transformation matrix between the source field mesh and the target field mesh is constructed by calling the radial basis function interpolation algorithm. The spatial interpolation mapping result of nodal forces and displacement increments on the global coupling interface is calculated using the spatial transformation matrix.

[0016] By adopting the above technical solution, solvers of different physical fields are allowed to perform calculations at their respective applicable time steps, thereby improving the overall solution efficiency. By setting synchronization barriers and fixed communication cycles, data exchange is carried out at designated nodes, while the radial basis function interpolation algorithm is used to perform data mapping between mismatched grids.

[0017] Preferably, in S4, the step of calculating the difference between the virtual work output by the source field and the virtual work received by the target field and accumulating them into a cumulative energy residual includes: Extract the source field nodal force vector output by the source field solver and the target field displacement increment vector output by the target field solver; The discrete virtual work transmitted from the source field physical domain to the global coupling interface is subtracted from the discrete virtual work received by the target field physical domain from the global coupling interface by the discrete inner product operation, and the discrete energy residual within the current preset macroscopic coupling communication cycle is calculated. The absolute value of the discrete energy residual is compared with the tolerance threshold. When the absolute value is less than the tolerance threshold, the discrete energy residual is algebraically summed with the historical residual value stored in the corresponding local energy residual accumulator, and the result is updated and written back to the local system to form the accumulated energy residual.

[0018] By adopting the above technical solution, energy exchange is calculated using virtual work, which has a clear physical meaning. By comparing it with the tolerance threshold, instantaneous numerical fluctuations below the threshold are avoided from being included in the cumulative error. This allows the accumulator to record the persistent energy deviation that leads to instability, providing quantitative data for subsequent corrections.

[0019] Furthermore, in S5, the step of determining the instability risk and triggering an interrupt signal based on the accumulated energy residual includes: In the established background monitoring thread, the absolute value of the accumulated energy residual is compared with the preset global tolerance limit in real time. When the absolute value exceeds the global tolerance limit, the running state transitions from the normal propulsion state to the alarm preparation state and the time window timer is started. The energy divergence rate is generated by calculating the average rate of change of the absolute value of the accumulated energy residual within a set sliding time window. The energy divergence rate is compared with the instability critical slope. When the energy divergence rate is greater than the instability critical slope, it is determined that the system has an instability risk, and the operating state transitions to the interrupt execution state, triggering the non-blocking type of interrupt signal.

[0020] By adopting the above technical solution, a dual criterion of absolute value and rate of change is established. The criterion based on the absolute value of the residual is used for total monitoring, and the criterion based on the energy divergence rate is used for trend analysis. This combination is used to distinguish between instantaneous error and continuously increasing error, so as to improve the accuracy of instability risk assessment.

[0021] Preferably, the instability critical slope is obtained through the following steps: The temporal energy fluctuation trajectories of different geometric configuration meshes running near the critical divergence boundary are extracted in advance, and the corresponding energy divergence rate features are extracted as training datasets. The state of whether or not a physical collapse is actually detected by the multiphysics simulation engine is defined as a system monitoring label; Based on the support vector machine classification algorithm, a classification hinge loss function is constructed between the predicted state and the true state. The decision boundary of the maximum margin hyperplane that distinguishes between stable and collapsed states is obtained by kernel function mapping and Lagrange multiplier optimization, which serves as the instability critical slope.

[0022] By adopting the above technical solution, this method determines the critical slope of instability through historical data and machine learning algorithms, replacing empirical manual settings. An SVM classifier is used to solve the classification boundary, and the generated instability criterion thus has data-driven generalization ability.

[0023] Further, in S6, the step of converting the accumulated energy residual into an additional nodal force vector by combining the spatial distance attenuation weight matrix includes: Calculate the Euclidean norm of the displacement increment vector of each node in the target field and construct the node displacement weight coefficient. Multiply the coefficient with the corresponding dimensionless weight value in the spatial distance attenuation weight matrix to construct the composite spatial distribution weight. The composite spatial distribution weight is used for reverse mapping and multiplied with the displacement increment direction vector corresponding to each node to generate the basic direction force vector. The cumulative energy residual is divided by the sum of the virtual work contributions of all basic direction force vectors to obtain the energy equivalent scaling factor, which is multiplied into the basic direction force vector to generate the initial error compensation load. Extract the discrete stiffness matrix and mass matrix corresponding to the target field control equation to establish a generalized eigenvalue solution equation, iteratively obtain a low-order rigid body mode vector set to construct an orthogonal projection operator, and perform spatial filtering on the initial error compensation load through the orthogonal projection operator to remove invalid energy components and generate the additional nodal force vector.

[0024] By adopting the above technical solution, this step constructs a method for converting energy and force. The composite spatial distribution weight makes the correction force mainly applied to nodes with large displacement and close to the boundary. The energy equivalent scaling makes the applied correction force correspond to the cumulative residual in terms of energy. The rigid body modes are filtered out by the orthogonal projection operator, so that the correction force mainly acts to eliminate deformation energy error, avoids the introduction of non-physical rigid body motion, and improves the stability of the correction process.

[0025] Preferably, in S6, the step of superimposing the external source terms into the target field solver includes: The generated additional node force vector is transformed into a dynamic source term function that decays according to an exponential function, wherein the time decay constant in the dynamic source term function is configured to be 0.5 to 2 times the preset macroscopic coupling communication period; The dynamic source term function is injected as an external load into the right-hand side of the non-homogeneous linear equation system of the target field solver to perform a source term reconstruction operation, and a synchronous release command is issued to unlock the synchronous barrier and restore the asynchronous time extrapolation process. After the action time window corresponding to the target field solver ends, the external load is removed to complete the energy smoothing correction, and the numerical data in the local energy residual accumulator is cleared to perform a cyclic iteration operation.

[0026] By adopting the above technical solution, this step applies the additional nodal force in the form of a function that decays smoothly over time, avoiding numerical shocks caused by step loads. The correction operation is coordinated with the time-progression loop of the solver. After the correction is completed, the accumulator is cleared so that the energy closed-loop correction system is executed cyclically during the simulation.

[0027] This invention provides a method for improving the accuracy of multiphysics coupled simulation solutions. It has the following beneficial effects: 1. This invention reduces cross-node memory access latency in large-scale parallel computing by directly associating the simulation region partitioning with the non-uniform memory access architecture of the underlying computing device. The method utilizes the affinity between the processor core and the physical memory segment to determine the data partitioning boundary, so that most of the data required by the computing task resides in its local memory in advance. This reduces communication overhead and waiting time in the solution process, and improves the overall computing efficiency and scalability of multiphysics coupled simulation.

[0028] 2. This invention establishes an instability risk prediction mechanism based on energy growth trends, changing the passive situation of traditional simulation methods that can only make judgments after the fact or rely on fixed thresholds. By calculating the rate of change of the cumulative energy residual in real time and comparing it with the instability critical slope obtained by training based on historical data, this method can identify the tendency of computational divergence in advance and trigger an interruption. This forward-looking judgment method avoids the continuous occurrence of invalid calculations and provides a valuable time window for timely intervention or preservation of effective intermediate results, thereby saving computing resources.

[0029] 3. This invention proposes a physically meaningful energy error correction scheme that transforms the accumulated energy residual into an equivalent additional nodal force. The spatial distribution of this force takes into account physical factors such as distance attenuation, and eliminates the influence of non-physical rigid body motion through orthogonal projection. The correction force is applied back to the model in the form of a smooth dynamic source term, ensuring the physical consistency of the correction process and avoiding the introduction of new numerical oscillations due to abrupt correction, thereby fundamentally improving the accuracy and reliability of simulation results. Attached Figure Description

[0030] Figure 1 This is a schematic diagram of the system architecture of the present invention; Figure 2 This is a schematic diagram of the method flow of the present invention; Figure 3 This is a schematic diagram illustrating the logic principle of local energy conservation domain partitioning and weight initialization in this invention; Figure 4 This is a schematic diagram of the timing of the dynamic calculation of the local tolerance threshold in this invention; Figure 5 This is a schematic diagram illustrating the asynchronous time integration and spatial mapping principle of the present invention; Figure 6 This is a schematic diagram of the discrete virtual work inner product calculation and residual local accumulation process of the present invention; Figure 7 This is a schematic diagram of the non-blocking interrupt determination flow topology of the present invention; Figure 8 This is a schematic diagram illustrating the principle of delayed orthogonal projection and source term reconstruction of the present invention; Figure 9This is a schematic diagram illustrating the trend of the extreme value of the interface energy residual as a function of the load fluctuation rate in this invention. Figure 10 This is a schematic diagram illustrating the trend of local node extreme displacement as a function of load fluctuation rate according to the present invention.

[0031] Among them, 10 is the processor; 20 is the memory; 30 is the communication bus; 101 is the hardware topology sensing and domain division module; 102 is the tolerance threshold calculation module; 103 is the asynchronous scheduling and space mapping module; 104 is the discrete virtual work calculation and state control module; and 105 is the delayed projection and source term reconstruction module. Detailed Implementation

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

[0033] Please see the appendix Figure 1 , Figure 1 This is a schematic diagram of the system architecture according to an embodiment of the present invention. The present invention provides a system for improving the accuracy of multiphysics coupling simulation, and provides a computing device hardware environment for supporting the operation of a method for improving the accuracy of multiphysics coupling simulation. The computing device may include: processor 10, memory 20 and communication bus 30.

[0034] The processor 10 establishes a low-level data transmission channel with the memory 20 through the communication bus 30. The computing device can run in a computing environment that supports multi-threaded concurrent execution. The memory 20 is configured with a non-uniform memory access architecture and is used to store multi-physics discrete mesh models and low-level data streams. The storage space of the memory 20 is divided into multiple independently addressable memory domains according to physical nodes. After the processor 10 reads the discrete mesh data in the memory 20, it extracts the mesh node coordinates and material properties for subsequent partial differential equation solving and physical quantity mapping processing.

[0035] The processor 10 is internally configured with various logic modules for executing the method of improving the accuracy of multiphysics coupling simulation. The logic modules may include: hardware topology perception and domain division module 101, tolerance threshold calculation module 102, asynchronous scheduling and spatial mapping module 103, discrete virtual work calculation and state control module 104, and delay projection and source term reconstruction module 105.

[0036] The hardware topology awareness and domain division module 101 is used to obtain the underlying memory physical distribution boundary, divide the global coupling interface into multiple local energy conservation domains, and calculate the spatial boundary attenuation weights corresponding to the generated grid nodes. The tolerance threshold calculation module 102 is used to extract the norm features of the governing equation matrix and the spatial gradient of the characteristic physical quantities in the local energy conservation domain, and to calculate and generate the tolerance threshold based on the matrix features and the physical quantity gradient. The asynchronous scheduling and spatial mapping module 103 is used to control the source field and target field solvers to perform deductions according to independent discrete time steps, and to perform spatial interpolation mapping of interface physical quantities at the set coupling communication points; The discrete virtual work calculation and state control module 104 is used to calculate the difference between the virtual work output by the interface source field and the virtual work received by the target field, accumulate the difference between the virtual work output by the interface source field and the virtual work received by the target field to the state register to form the accumulated energy residual, and trigger an interrupt signal based on the comparison result of the accumulated energy residual with the global tolerance extreme value and the instability critical slope. The delayed projection and source term reconstruction module 105 is used to perform orthogonal projection by combining the spatial boundary attenuation weight and the accumulated energy residual when an interruption signal is detected, reconstruct the additional node force vector and superimpose it into the external source term injected into the target field.

[0037] The method for improving the accuracy of multiphysics coupled simulation solutions described below can be referred to in correspondence with the system for improving the accuracy of multiphysics coupled simulation solutions described above.

[0038] The method in this embodiment can be used to execute the above system embodiment, and its principle and technical effect are similar, so it will not be described again here.

[0039] Please see the appendix Figure 2 , Figure 2 This is a schematic flowchart of a method according to an embodiment of the present invention. The present invention provides a method for improving the accuracy of multiphysics coupling simulation solutions, comprising the following steps: S10, the hardware topology awareness and domain division module 101 extracts the physical node boundary of the non-uniform memory access architecture in the memory 20, divides the global coupling interface into local energy conservation domains residing in independent memory addresses based on the physical node boundary, and calculates and generates the spatial distance attenuation weight matrix on the cross-domain boundary. S20, the tolerance threshold calculation module 102 extracts the numerical stiffness feature norm and the spatial gradient extremum of the feature physical quantity corresponding to each local energy conservation domain, and calls the numerical stiffness feature norm and the spatial gradient extremum of the feature physical quantity to calculate and generate the tolerance threshold of local numerical stability. S30, the asynchronous scheduling and spatial mapping module 103 sends instructions to the source field solver and the target field solver, controlling the source field solver and the target field solver to perform asynchronous time advancement according to their respective discrete time steps, and to perform spatial interpolation mapping of nodal forces and displacement increments of the local energy conservation domain interface in parallel when the time axis is extrapolated to the corresponding coupled communication point. S40, Discrete Virtual Work Calculation and State Control Module 104 extracts the nodal force and displacement increments of the source field and the target field, calculates the difference between the virtual work output by the source field and the virtual work received by the target field in the current communication cycle through discrete inner product operation, and accumulates the difference between the virtual work output by the source field and the virtual work received by the target field in the current communication cycle into the state register corresponding to the local energy conservation domain to form the accumulated energy residual. S50, the discrete virtual work calculation and state control module 104 determines the system operating state based on the accumulated energy residual and its rate of change over time in the state register, combined with the preset global tolerance extreme value and instability critical slope criterion; when the determination result shows that the accumulated energy residual has the risk of instability, a non-blocking interrupt signal is triggered and the current integration process continues to run to the next coupling communication point; otherwise, when the determination result shows that the system is in a stable state, no additional operation is required and the solver is instructed to continue to execute the next time integration step; S60, the delayed projection and source term reconstruction module 105 loads the rigid body modal orthogonal projection operator and spatial distance attenuation weight matrix for the local energy conservation domain that triggers the interrupt signal when assembling the initial boundary conditions for the next time step. It converts the accumulated energy residual into an additional nodal force vector, superimposes the additional nodal force vector onto the basic source term load vector of the target field solver, and after the superposition is completed, notifies the discrete virtual work calculation and state control module 104 to clear the numerical data in the state register to perform the loop iteration operation.

[0040] See attached document Figure 3 , Figure 3 This is a schematic diagram of the logic principle of local energy conservation domain partitioning and weight initialization according to an embodiment of the present invention. In this embodiment, the processor 10 calls the underlying interface instruction to drive the hardware topology sensing and domain division module 101 to perform the exploration of the underlying physical architecture and the reorganization and decomposition of the spatial mesh. The hardware topology sensing and domain division module 101 extracts the physical node boundaries of the non-uniform memory access architecture in the memory 20. The hardware topology sensing and domain division module 101 divides the global coupling interface into local energy conservation domains residing in independent memory addresses according to the physical node boundaries. The hardware topology sensing and domain division module 101 calculates and generates the spatial distance attenuation weight matrix on the cross-domain boundary.

[0041] In the specific implementation process, the local energy conservation domain partitioning and weight initialization process based on hardware topology reads the non-uniform memory access architecture configuration table through the application programming interface of the operating system kernel installed on the computing device. The hardware topology perception and domain division module 101 parses the non-uniform memory access architecture configuration table to extract the affinity distance matrix between the processing core inside the processor 10 and different physical addressing segments of the memory 20. The hardware topology perception and domain division module 101 determines the physical node boundary at the bottom layer of the computing device based on the affinity distance matrix. The hardware topology perception and domain division module 101 establishes a memory access latency mapping map in a multi-threaded concurrent environment by integrating the physical node boundary information at the bottom layer.

[0042] As a preferred approach, the hardware topology sensing and domaining module 101 loads the global coupling interface mesh node data from the multiphysics discrete mesh model. Based on the mesh node coordinate topological relationship and memory access delay mapping, the hardware topology sensing and domaining module 101 calls a multi-level graph partitioning algorithm to process the mesh data. The layer coarsening and boundary segmentation operations of the multi-level graph partitioning algorithm can be completed using existing multi-level graph partitioning algorithm resource libraries. The specific implementation process of the multi-level graph partitioning algorithm will not be described in detail here. The hardware topology sensing and domaining module 101 uses the multi-level graph partitioning algorithm to divide the global coupling interface into multiple non-overlapping local energy conservation domains that are spliced ​​and cover the global coupling interface. The hardware topology sensing and domaining module 101 allocates and locks the mesh node coordinates and material property data contained in each local energy conservation domain to the corresponding independent memory address partitioned by the physical node boundary.

[0043] The hardware topology sensing and domain delineation module 101 establishes a boundary distance decay kernel function mathematical model for each set of grid nodes contained within a local energy conservation domain. The specific methods for creating and training this model are provided below. The hardware topology sensing and domain delineation module 101 collects historical perturbation propagation datasets of cross-domain interfaces under different grid sizes as training samples. The hardware topology sensing and domain delineation module 101 extracts the optimal reference weight distribution that minimizes the local energy residual at the interface and defines it as a supervision label. The hardware topology sensing and domain delineation module 101 constructs a loss function based on the mean square error between the predicted weight distribution and the supervision label. The hardware topology sensing and domain delineation module 101 calls a nonlinear least squares optimization algorithm to iteratively update the spatial decay feature length constant within the model until the loss function converges (e.g., the loss function value is less than 1E-6), thereby solidifying the decay control parameters of the boundary distance decay kernel function mathematical model. In this embodiment, the typical range of the spatial decay feature length constant λ is 2 to 10 times the average grid size within the local energy conservation domain.

[0044] After parameter training is completed, the boundary distance decay kernel function mathematical model exhibits an adjustment mechanism that decreases exponentially with spatial span to smooth the jumps in physical field values ​​at the boundary. The hardware topology sensing and domain delineation module 101 traverses any discrete grid node within the local energy conservation domain and calculates the shortest Euclidean distance scalar from any discrete grid node to the corresponding cross-domain boundary. The hardware topology sensing and domain delineation module 101 substitutes the shortest Euclidean distance scalar into the calculation to generate the corresponding dimensionless weight value. The boundary distance decay kernel function mathematical model calculated and called by the hardware topology sensing and domain delineation module 101 is expressed as follows: ; In the formula, For the first in the local energy conservation domain Dimensionless weight values ​​of discrete grid nodes; This refers to the index number of the discrete grid nodes within the local energy conservation domain; For the natural constant The natural exponential function operator with base 0; A scalar operator for calculating the shortest Euclidean distance between two spatial geometric objects; For the first The three-dimensional spatial coordinate vector of the discrete grid node in the global coordinate system; It is the cross-domain boundary where a local energy-conserving domain and an adjacent energy-conserving domain are in close contact; The length constant of the spatial decay feature is obtained through data-driven training.

[0045] After completing the calculation of the dimensionless weight values ​​of all discrete grid nodes in the current local energy conservation domain, the hardware topology sensing and domain delineation module 101 constructs a matrix-form data structure. The hardware topology sensing and domain delineation module 101 constructs a diagonal matrix-form spatial distance attenuation weight matrix using all dimensionless weight values ​​as the main diagonal elements of the matrix. The hardware topology sensing and domain delineation module 101 persists the spatial distance attenuation weight matrix to the independent memory address corresponding to the local energy conservation domain so that the delayed projection and source term reconstruction module 105 can initiate a call request when reconstructing the force vector of the additional node.

[0046] See attached document Figure 4 , Figure 4This is a schematic diagram of the timing of dynamic calculation of local tolerance threshold according to an embodiment of the present invention. In this embodiment, the processor 10 schedules the tolerance threshold calculation module 102 to perform real-time evaluation of local calculation state and generation of error judgment benchmark. The tolerance threshold calculation module 102 avoids using globally uniform fixed parameters during the physical field time integration derivation. The tolerance threshold calculation module 102 extracts the norm features of the control equation matrix and the spatial gradient of the characteristic physical quantities in the local energy conservation domain. The tolerance threshold calculation module 102 dynamically calculates and generates the tolerance threshold for a single independent addressable memory region based on the matrix features and the physical quantity gradient.

[0047] In the specific implementation process, the dynamic calculation mechanism of the local tolerance threshold is initiated by the tolerance threshold calculation module 102 extracting the system Jacobian matrix or stiffness matrix generated after discretizing the corresponding local energy conservation domain control equation within the current time integration step. From the perspective of physical principles, the system Jacobian matrix or stiffness matrix directly reflects the inherent physical characteristics of the local mesh structure resisting numerical disturbances. The tolerance threshold calculation module 102 calls the matrix norm calculation operator to perform Frobenius norm calculation on the extracted system Jacobian matrix or stiffness matrix. The process of taking the square root of the sum of squares of matrix elements in the Frobenius norm operation can be completed using existing numerical linear algebra techniques. The tolerance threshold calculation module 102 generates a numerical stiffness characteristic norm that can characterize the overall numerical stiffness of the local energy conservation domain through norm calculation.

[0048] As a preferred approach, the tolerance threshold calculation module 102 extracts the numerical stiffness feature norm and then locks the main control feature physical quantity in the multiphysics simulation variables. The main control feature physical quantity is determined according to the specific coupled simulation scenario and configured as flow field pressure or solid displacement variables in the fluid-structure interaction scenario. The tolerance threshold calculation module 102 reads the values ​​of the main control feature physical quantity of all discrete grid nodes inside the local energy conservation domain and calculates the physical quantity gradient vector of the discrete grid nodes in three-dimensional space based on the finite volume or finite element shape function interpolation algorithm. The tolerance threshold calculation module 102 calculates the second norm of the physical quantity gradient vector of all discrete grid nodes and performs traversal comparison. The tolerance threshold calculation module 102 selects the gradient second norm with the largest value inside the local energy conservation domain as the spatial gradient extremum of the feature physical quantity to capture the local stress concentration phenomenon caused by mesh incoordination or material mutation in the interface region.

[0049] After acquiring local feature parameters, the tolerance threshold calculation module 102 fuses the two to construct a dynamic tolerance threshold calculation mathematical model. This mathematical model also relies on data-driven production and training. The tolerance threshold calculation module 102 extracts the stiffness norm and gradient extrema from the historical multiphysics coupling deduction process as sample inputs. The tolerance threshold calculation module 102 records whether the simulation system eventually experiences numerical collapse as a binary discrete supervision label. Based on the input and label, the tolerance threshold calculation module 102 constructs a hinge loss function that includes a regularization penalty term and aims to maximize the classification boundary between stable and divergent samples. The tolerance threshold calculation module 102 uses an adaptive moment estimation gradient descent algorithm to perform multiple rounds of parameter optimization on the empirical scaling constant inside the model along the opposite direction of the loss function gradient until the loss function converges and the optimal calculation coefficients of this tolerance judgment model are locked. In this embodiment, the typical range of the empirical scaling constant α is 1E-4 to 1E-2.

[0050] The tolerance threshold calculation module 102 multiplies the matrix stiffness feature with the spatial gradient extremum through multiplication and substitutes it into a trained and fixed scaling constant to unify the dimensions, thereby calculating and generating a tolerance threshold for local numerical stability. The tolerance threshold calculation module 102 writes the tolerance threshold into a dedicated memory register of the corresponding local energy conservation domain and reactivates the calculation operation every time a macroscopic coupling communication cycle is advanced. The computing device adaptively adjusts the error judgment boundary through this mechanism. The mathematical model called by the tolerance threshold calculation module 102 is expressed as the following formula: ; In the formula, In order to target the The tolerance threshold is dynamically calculated and generated within a local energy conservation domain; The index number of the local energy-conserving domain generated by the global partitioning; These are empirical scaling constants that have been solidified through model parameter optimization; Calculate operators for the Frobenius norm of a matrix; For the first The system Jacobian matrix or stiffness matrix generated after discretizing the system according to the governing equations of a local energy conservation domain; To define a domain in a specified continuous space Operators for extracting extrema from within the core; A vector representing the coordinates of any spatial location within the domain of a continuous space; For the first The continuous spatial domain corresponding to a local energy conservation domain; Calculate operators for the 2-norm (Euclidean norm) of spatial vectors; It is a three-dimensional gradient operator; For spatial position coordinate vectors The main control feature physical quantity extracted from the site.

[0051] See attached document Figure 5 , Figure 5 This is a schematic diagram of the asynchronous time integration and spatial mapping principle according to an embodiment of the present invention. In this embodiment, the processor 10 calls the asynchronous scheduling and spatial mapping module 103 to issue instructions to the source field solver and the target field solver. The asynchronous scheduling and spatial mapping module 103 controls the source field solver and the target field solver to perform asynchronous time advancement according to independent discrete time steps. As a preferred method, the independent discrete time step is calculated and generated according to the Courant-Friedrich-Lévy stability condition. In this embodiment, the Courant number is preferably in the range of 0.1 to 0.8.

[0052] The calculation process extracts the minimum grid size and medium wave velocity of the source field and target field, respectively. Under the premise of ensuring that the maximum step size allowed by the Courant-Friedrich-Levy stability condition is not exceeded, the asynchronous scheduling and spatial mapping module 103 configures the time step size of the source field and the time step size of the target field to be an integer fraction of the macroscopic coupling communication period. In this embodiment, the integer is preferably a value between 10 and 1000 to ensure that the local advance clock can be accurately aligned with the coupling communication point on the global time axis. The macroscopic coupling communication period is a global synchronization interval that is preset according to the characteristic time scale of the coupled physical phenomenon.

[0053] The asynchronous scheduling and spatial mapping module 103 continuously monitors the local advance clocks of the source field solver and the target field solver within the time integration cycle. The source field solver and the target field solver independently solve the control equations and update the physical quantity states in their respective sub-time steps. When the local advance clock of the source field solver or the target field solver has not reached the macroscopic coupling communication cycle, the asynchronous scheduling and spatial mapping module 103 instructs the corresponding solver to continue executing the internal time step advance. When the local advance clocks of the source field solver and the target field solver have both evolved to the corresponding coupling communication point, the asynchronous scheduling and spatial mapping module 103 triggers a synchronization barrier instruction. The asynchronous scheduling and spatial mapping module 103 suspends the time evolution process of the source field solver and the target field solver through the synchronization barrier instruction and enters the spatial data exchange stage.

[0054] During the spatial data exchange phase, the asynchronous scheduling and spatial mapping module 103 performs physical quantity transfer operations on the non-coordinated grid nodes on both sides of the local energy conservation domain interface. The asynchronous scheduling and spatial mapping module 103 extracts the boundary node force vector output by the source field solver and the displacement increment vector output by the target field solver. The spatial mapping operation based on the non-matching grid interface is explained from the physical principle level as smoothly interpolating the physical quantity field discrete on the source field grid to the target field grid nodes through basis functions. The asynchronous scheduling and spatial mapping module 103 calls the radial basis function interpolation algorithm to construct the spatial transformation matrix between the source field grid and the target field grid. As a preferred method, the kernel function of the radial basis function interpolation algorithm can be a multi-quadrant function or a Gaussian function. The solution process of the linear equation system of the radial basis function interpolation algorithm can be completed based on the existing scattered data spatial interpolation technology.

[0055] The asynchronous scheduling and spatial mapping module 103 calculates the spatial interpolation mapping results of nodal forces and displacement increments on the non-coordinated mesh interface through the spatial transformation matrix. The spatial interpolation mapping mathematical model called by the asynchronous scheduling and spatial mapping module 103 is expressed as follows: ; ; In the formula, This is the target field node force vector mapped to the target field grid nodes through spatial interpolation; This is the force transformation interpolation matrix projected from the source field grid to the target field grid; This is the source field nodal force vector output by the source field solver at the current coupled communication point; This is the source field displacement increment vector mapped to the source field grid nodes via spatial interpolation; This is the displacement transformation interpolation matrix projected from the target field grid to the source field grid; This is the target field displacement increment vector output by the target field solver at the current coupled communication point.

[0056] See attached document Figure 6 , Figure 6This is a schematic diagram of the discrete virtual work inner product calculation and residual local accumulation process according to an embodiment of the present invention. In this embodiment, the processor 10 calls the discrete virtual work calculation and state control module 104 to perform the energy conservation verification operation of the multi-physics coupling interface. The discrete virtual work calculation and state control module 104 receives the source field node force vector output by the source field solver and the target field displacement increment vector output by the target field solver from the output buffer of the asynchronous scheduling and space mapping module 103. At the same time, the discrete virtual work calculation and state control module 104 retrieves the force transformation interpolation matrix and displacement transformation interpolation matrix generated by the asynchronous scheduling and space mapping module 103. Based on the received physical quantity vector and interpolation matrix, the discrete virtual work calculation and state control module 104 calculates the discrete energy residual introduced by asynchronous propulsion and space mapping in the current macroscopic coupling communication cycle.

[0057] In the specific implementation process, the physical principle of the discrete virtual work inner product calculation operation is that, under the ideal energy conservation condition, the work transferred from the source field perspective (the work done by the source field force on the target field displacement after mapping) should be equal to the work received from the target field perspective (the work done by the source field force on the target field displacement after mapping). Any deviation introduced by the numerical algorithm will be reflected as the difference between the two being non-zero. The discrete virtual work calculation and state control module 104 obtains the discrete energy residual that can accurately characterize the degree of energy non-conservation by subtracting the discrete virtual work transferred from the source field physical domain to the coupling interface from the discrete virtual work received by the target field physical domain from the coupling interface.

[0058] As a preferred approach, after calculating the discrete energy residual, the discrete virtual work calculation and state control module 104 reads the dynamic tolerance threshold generated by the tolerance threshold calculation module 102 from the dedicated memory register of the corresponding local energy conservation domain. The discrete virtual work calculation and state control module 104 compares the absolute value of the discrete energy residual with the dynamic tolerance threshold. When the absolute value of the discrete energy residual is less than the dynamic tolerance threshold, the energy imbalance generated in the current coupling communication cycle is considered acceptable. At this time, the discrete virtual work calculation and state control module 104 sends the value of the discrete energy residual to the local energy residual accumulator for storage and accumulation.

[0059] As a protection mechanism for extreme cases, when the absolute value of the discrete energy residual is much greater than the dynamic tolerance threshold (e.g., more than 5 times the dynamic tolerance threshold), an immediate coupled iterative correction process can be selectively triggered. For example, by adjusting the relaxation factor of the boundary conditions or requesting the asynchronous scheduling and space mapping module 103 to perform a rollback and recalculation of a sub-time step, to prevent single-step calculation divergence. In this embodiment, the relaxation factor is preferably in the range of 0.5 to 0.9. Even if the absolute value of the discrete energy residual is slightly greater than or equal to the dynamic tolerance threshold, the discrete virtual work calculation and state control module 104 still performs the residual accumulation operation and hands over the energy correction decision to the subsequent state machine.

[0060] The mathematical model for calculating the discrete energy residual executed by the discrete virtual work calculation and state control module 104 is expressed as follows: ; In the formula, Characterization in the Discrete energy residuals calculated at the interface of a local energy conservation domain during the current coupled communication cycle; This is the source field nodal force vector output by the source field solver at the current coupled communication point; An algebraic transpose operator performed on a vector or matrix; This is the displacement transformation interpolation matrix projected from the target field grid to the source field grid; This is the target field displacement increment vector output by the target field solver at the current coupled communication point; is the force transformation interpolation matrix projected from the source field grid to the target field grid.

[0061] After residual verification or iterative correction, the discrete virtual work calculation and state control module 104 writes the finally determined discrete energy residual into the local energy residual accumulator bound to the current local energy conservation domain. The local energy residual accumulator is physically implemented as a dedicated memory address. The accumulation operation performed by the discrete virtual work calculation and state control module 104 is to perform algebraic summation between the currently calculated discrete energy residual and the historical residual value stored in the local energy residual accumulator, and write the updated accumulated energy residual back to the local energy residual accumulator. The accumulated energy residual will be used to construct the error compensation load in subsequent calculation cycles.

[0062] See attached document Figure 7 , Figure 7This is a schematic diagram of the non-blocking interrupt judgment flow topology according to an embodiment of the present invention. In this embodiment, the processor 10 calls the discrete virtual work calculation and state control module 104 to establish an independent background monitoring thread within the multiphysics parallel computing framework. The discrete virtual work calculation and state control module 104 continuously monitors the accumulated energy residual value in the local energy residual accumulator in the background monitoring thread. The discrete virtual work calculation and state control module 104 constructs a non-blocking judgment logic with multiple discrete operating states based on the finite state machine theory. The discrete virtual work calculation and state control module 104 evaluates whether the accumulated energy residual exceeds the set system safety limit without suspending the main operation processes of the origin field solver and the target field solver.

[0063] In the specific implementation process, the non-blocking interrupt judgment flow logic relies on the state transition mechanism of the finite state machine to realize the lock-free management of complex asynchronous events. The discrete virtual power calculation and state control module 104 is pre-configured with normal progress state, alarm preparation state and interrupt execution state. The normal progress state indicates that the current accumulated energy residual is within the tolerable calculation noise range of the multiphysics simulation system. The discrete virtual power calculation and state control module 104 compares the absolute value of the accumulated energy residual in the local energy residual accumulator with the preset global tolerance extreme value in real time.

[0064] In this embodiment, the global tolerance limit is preferably configured to be 5 to 10 times the average value of the dynamic tolerance thresholds corresponding to all local energy conservation domains. When the absolute value of the accumulated energy residual exceeds the global tolerance limit, the operating state of the discrete virtual work calculation and state control module 104 transitions from the normal propulsion state to the alarm preparation state. The discrete virtual work calculation and state control module 104 in the alarm preparation state starts a time window timer and records the number of continuous communication steps exceeding the global tolerance limit. In this embodiment, the number of continuous communication steps is set to 3 to 10 consecutive communication steps. When the absolute value of the accumulated energy residual falls back below the global tolerance limit, the discrete virtual work calculation and state control module 104 returns from the alarm preparation state to the normal propulsion state and clears the time window timer.

[0065] When the discrete virtual work calculation and state control module 104 is in the alarm preparation state, it calls the mathematical model for calculating the energy divergence rate and combines it with the critical slope of instability to determine the physical instability phenomenon. The method of making and training the key parameters of the determination model is disclosed here. The discrete virtual work calculation and state control module 104 pre-extracts the time-series energy fluctuation trajectory of different geometric configuration grids running near the critical divergence boundary and calculates and extracts the corresponding energy divergence rate features as training datasets.

[0066] The Discrete Virtual Work Calculation and State Control Module 104 defines the state of physical collapse actually detected by the simulation engine as the system supervision label. The Discrete Virtual Work Calculation and State Control Module 104 constructs a classification hinge loss function between the predicted state and the real state based on the support vector machine classification algorithm. The Discrete Virtual Work Calculation and State Control Module 104 solves the maximum margin hyperplane that distinguishes between stable and collapsed states by kernel function mapping and Lagrange multiplier optimization. Based on the data distribution, the decision boundary of the hyperplane is extracted as the instability critical slope constant for the final state assessment.

[0067] The Discrete Virtual Work Calculation and State Control Module 104 calculates the average rate of change of the absolute value of the accumulated energy residual within the sliding time window based on physical principles. The mathematical model for calculating the energy divergence rate called by the Discrete Virtual Work Calculation and State Control Module 104 is expressed as follows: ; In the formula, To monitor the calculated energy dissipation rate within the sliding time window; In this embodiment, the total number of macroscopic coupling communication cycles covered by the set sliding time window is... The typical value range is 5 to 50; It is a summation operator for discrete sequences; This is the global index number of the current macroscopic coupling communication cycle in the physical field deduction process; This is a cyclic index variable for the macroscopic coupling communication cycle within the sliding time window; Operators for extracting the absolute value of numerical values; In the first Local cumulative energy residuals extracted from macroscopic coupled communication points; In the first Local cumulative energy residuals extracted from macroscopic coupled communication points; This is the global time extrapolation interval between two adjacent macroscopically coupled communication cycles.

[0068] After calculating the energy divergence rate, the discrete virtual work calculation and state control module 104 compares the energy divergence rate with the system's preset instability critical slope. When the energy divergence rate is greater than the instability critical slope, the discrete virtual work calculation and state control module 104 transitions from the alarm preparation state to the interrupt execution state. The instability critical slope is calibrated based on historical simulation experience data and configured as a positive real number in this embodiment. In this embodiment, the typical value range of the instability critical slope is usually in the range of 1.0 to 10.0 after normalization based on the dimensions of the energy residual. The discrete virtual work calculation and state control module 104, which is in the interrupt execution state, sends an asynchronous interrupt signal containing the accumulated energy residual value to the delay projection and source term reconstruction module 105, and simultaneously issues a barrier interception command to the asynchronous scheduling and space mapping module 103.

[0069] After receiving the interception command, the asynchronous scheduling and spatial mapping module 103 forces the source field solver and the target field solver to enter the synchronization barrier at the next nearest macroscopic coupling communication point to cooperate with the subsequent energy correction operation of the delay projection and source term reconstruction module 105.

[0070] After the delay projection and source term reconstruction module 105 completes the energy correction operation, the discrete virtual work calculation and state control module 104 clears the local energy residual accumulator and resets the running state to the normal propagation state. The discrete virtual work calculation and state control module 104 restores the non-blocking background monitoring of the multiphysics time integration derivation process through the state reset mechanism. The allocation of the state register of the finite state machine and the creation of the background independent thread can be completed based on the existing concurrent programming and microprocessor architecture technology. The above-mentioned non-blocking interrupt judgment and flow logic gives the multiphysics parallel computing framework the ability to adaptively recover when dealing with local numerical disturbances and avoids unnecessary suspension of global computing resources.

[0071] See attached document Figure 8 , Figure 8 This is a schematic diagram of the delayed orthogonal projection and source term reconstruction principle according to an embodiment of the present invention. In this embodiment, the processor 10 calls the delayed projection and source term reconstruction module 105 to perform an energy correction operation in response to the asynchronous interrupt signal issued by the discrete virtual work calculation and state control module 104. The delayed projection and source term reconstruction module 105 extracts the accumulated energy residual transmitted by the interrupt signal and the target field displacement increment vector obtained from the target field.

[0072] Based on the equivalence principle of the right-hand side of the multiphysics discrete control equation, the scalar form of the accumulated energy residual is distributed to each node of the non-coordinated mesh interface. The delayed projection and source term reconstruction module 105 calculates the Euclidean norm of the displacement increment vector of each node in the target field and constructs the proportion coefficient of the node displacement relative to the sum of the norms of the displacement increments of all target field nodes in the local energy conservation domain. At the same time, it sends a call request to the corresponding independent memory address to read the spatial distance attenuation weight matrix generated by the hardware topology perception and domain division module 101. The delayed projection and source term reconstruction module 105 multiplies the above proportion coefficient with the corresponding dimensionless weight value in the spatial distance attenuation weight matrix to construct the composite spatial distribution weight.

[0073] The delayed projection and source term reconstruction module 105 uses the composite spatial distribution weight to perform reverse mapping and multiplies the composite spatial distribution weight with the displacement increment direction vector corresponding to each node to generate the basic direction force vector. Then, it calculates the test virtual work done by the entire basic direction force vector on the corresponding target field displacement increment, and uses the quotient of the accumulated energy residual and the test virtual work as the energy equivalent scaling factor. The energy equivalent scaling factor is multiplied into the basic direction force vector to ensure that the work done by the generated force vector strictly matches the energy error, thereby generating the initial error compensation load in the form of local node force vector.

[0074] The delayed projection and source term reconstruction module 105 performs spatial filtering on the initial error compensation load through the mathematical model constructed by the delayed orthogonal projection mechanism. The production of this orthogonal projection operator model requires the pre-execution of modal feature extraction and data training processes. The delayed projection and source term reconstruction module 105 extracts the discrete stiffness matrix and mass matrix corresponding to the target field control equation. The delayed projection and source term reconstruction module 105 does not rely on external dynamic response data, but directly establishes generalized eigenvalue solution equations based on the inherent characteristics of the system itself. The delayed projection and source term reconstruction module 105 uses the subspace iteration method or Lanczos algorithm to repeatedly update the approximation vector of the stiffness and mass feature matrix in the implicit space to obtain the low-order rigid body modal vector set, thereby finally assembling the operator model that conforms to the physical orthogonality characteristics.

[0075] The orthogonal projection operation mathematically removes the invalid energy components parallel to the low-order rigid body modes from the initial error compensation load. The mathematical model of delayed orthogonal projection called by the delayed projection and source term reconstruction module 105 is expressed as follows: ; In the formula, Characterizes the effective compensation load vector generated after orthogonal projection filtering; The identity matrix that matches the dimension of the initial error compensation load vector; The low-order rigid body mode matrix is ​​characterized by orthogonal normalization and the column vectors of the low-order rigid body mode matrix are linearly independent in the spatial domain. The transpose of the low-order rigid body mode matrix; The initial error compensation load vector is characterized by the back mapping of the accumulated energy residual and the target field displacement increment vector through spatial distribution weights.

[0076] After completing the orthogonal projection operation, the delayed projection and source term reconstruction module 105 executes the external source term compensation and reconstruction logic. The delayed projection and source term reconstruction module 105 converts the effective compensation load vector into a dynamic source term function that decays according to an exponential function. The time decay constant in the dynamic source term function is preferably configured to be 0.5 to 2 times the macroscopic coupling communication period in order to achieve a smooth release of the correction energy.

[0077] The delayed projection and source term reconstruction module 105 injects the dynamic source term function as an external load into the right-hand side of the non-homogeneous linear equation system of the target field solver to perform source term reconstruction. Then, it sends a synchronization release command to the asynchronous scheduling and spatial mapping module 103 to unlock the synchronization barrier and restore the asynchronous time extrapolation process of the multiphysics simulation system. The target field solver continuously calculates the physical state including the energy correction effect in the local time step corresponding to the subsequent action time window. After the action time window ends, the delayed projection and source term reconstruction module 105 cancels the external load to complete the energy smoothing correction. The extraction of the low-order rigid body mode matrix can be completed using existing finite element modal analysis techniques.

[0078] To aid in understanding the technical solution of this invention, an application example of aerodynamic-thermal structure coupling simulation of an aero-engine turbine blade is provided below.

[0079] This embodiment applies the accuracy improvement system of multiphysics coupling simulation to bidirectional extrapolation calculations involving gas flow field and metal blade solid field. The computing device motherboard is configured with a non-unified memory access architecture. The system divides the main control memory into 4 independently addressable physical nodes according to the underlying hardware manual. The macroscopic coupling communication cycle is set to 100 microseconds based on the time characteristic scale of fluid mechanics and solid mechanics physical fields.

[0080] In the initial stage of system operation, the tolerance threshold calculation module 102 extracts the pressure spatial gradient of the fluid grid nodes and the stiffness matrix of the solid grid structure and generates the corresponding tolerance threshold of the leading edge local energy conservation domain as 4.5E-3 Joules. The asynchronous scheduling and spatial mapping module 103 issues independent derivation commands to control the flow field solver to run at a time step of 0.1 microseconds and the solid field solver to run at a time step of 5 microseconds.

[0081] When the simulation reaches the 850th macroscopic coupling communication cycle, the discrete virtual work calculation and state control module 104 calculates the difference between the virtual work output by the source field and the virtual work received by the target field in the leading edge local energy conservation domain through discrete inner product operation. The difference is 1.2E-2 joules. This value exceeds the local tolerance threshold of 4.5E-3 joules and accumulates to 5.8E-2 joules in the following 4 communication cycles. The discrete virtual work calculation and state control module 104 calculates the energy divergence rate to be 1.15, which is greater than the set instability critical slope of 1.0. Based on this, the system triggers a non-blocking type interrupt signal.

[0082] When the system triggers the interrupt execution state, the delayed projection and source term reconstruction module 105 retrieves the spatial distance attenuation weight matrix of the corresponding region and reverse maps the cumulative energy residual of 5.8E-2 joules to the basic direction force vector. The module performs modal orthogonal projection operation to filter out rigid body translation components, reconstructs additional nodal force vectors and superimposes them into the external source term of the solid field solver in the form of an exponential function with a time decay constant of 150 microseconds. After the reconstruction is executed, the local energy residual value of the interface decreases and the state progresses.

[0083] To verify the technical performance of the above scheme under different aerodynamic load fluctuation conditions in high-performance computing nodes, a control group using a fixed relaxation factor synchronization strategy and an experimental group using the accuracy improvement strategy of this invention were introduced.

[0084] Table 1: Performance Comparison Test Data of Traditional Coupled Solution Strategy and Accuracy Improvement Strategy According to Table 1 and Figure 9 and Figure 10 It can be seen that the interface energy residual value measured by the control group using the fixed coupling solution strategy under various aerodynamic load fluctuation conditions increases monotonically from 5.781 joules to 7.853 joules as the aerodynamic load fluctuation rate increases. The energy residual exceeding the tolerance causes the solver to diverge and restart 8 to 16 times per 10,000 steps. The local node extreme displacement of the control group is distributed in the range of 18.42 mm to 28.24 mm and exceeds the set limit. Conventional time integration and spatial interpolation operations generate energy injection at the non-coordinated mesh interface.

[0085] After introducing the accuracy improvement system and method of this invention, the extreme value of the interface energy residual in the experimental group decreased and remained in the range of 0.035 joules to 0.048 joules. The extreme displacement of local nodes fell back to the range of 1.52 mm to 1.76 mm and fell within the safety constraints generated by the tolerance threshold calculation module 102. The number of divergence restarts of the system solver decreased to 0. After performing the additional node force vector reconstruction operation, the average single-step synchronization time of the experimental group was in the range of 2.41 seconds to 2.58 seconds, which was higher than the range of 2.15 seconds to 2.24 seconds of the control group. The delay projection and source term reconstruction module 105 occupied additional calculation instruction cycles when calculating the orthogonal projection operator and reconstructing the external source term and complied with the global blocking constraint boundary.

[0086] The accuracy enhancement mechanism in this scheme has the ability to coordinate intervention of multiple physics fields. The spatial distance attenuation weight matrix established by the hardware topology perception and domain division module 101 maps the underlying physical memory distribution attributes. The tolerance threshold calculation module 102, together with the discrete virtual work calculation and state control module 104, quantifies the local energy residual and divergence slope. The system filters out translation components within the set eigenvalue boundary by executing non-blocking interrupt signal issuance and delayed orthogonal projection. The system avoids the accumulation of discrete errors in the same spatial location and realizes physical state intervention in the process of solving the multi-physics field coupled equations.

[0087] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for improving the accuracy of multiphysics coupled simulation solutions, wherein the method is configured to be executed in a computing device, characterized in that, The computing device includes a processor, a memory, and a communication bus. The processor establishes a low-level data transmission channel with the memory through the communication bus. The method includes the following steps: S1, the global coupling interface is divided into a local energy conservation domain according to the hardware topology boundary of the computing device, and a spatial distance attenuation weight matrix is ​​generated. S2, extract the numerical stiffness feature norm and the spatial gradient extrema of the feature physical quantity of the local energy conservation domain to calculate and generate the tolerance threshold. S3 controls the source field solver and the target field solver to advance asynchronously and performs spatial interpolation mapping at the coupling communication point; S4, calculate the difference between the virtual work output by the source field and the virtual work received by the target field and accumulate it as the cumulative energy residual; S5, determine the instability risk based on the accumulated energy residual and trigger an interrupt signal; S6, combining the spatial distance attenuation weight matrix, the accumulated energy residual is transformed into an additional nodal force vector, which is then superimposed on the external source term of the target field solver.

2. The method for improving the accuracy of multiphysics coupled simulation solutions according to claim 1, characterized in that, In S1, the step of dividing the global coupling interface into a local energy-conserving domain according to the underlying hardware topology boundary of the computing device includes: Extract non-uniform memory access architecture configuration information, parse and extract the affinity distance matrix between the processor's internal processing core and different physical addressing segments of the memory, and determine the physical node boundary at the bottom layer of the computing device based on the affinity distance matrix; The multi-level graph partitioning algorithm is called to divide the global coupling interface into multiple non-overlapping and spliced ​​local energy conservation domains, and the grid node coordinates and material property data contained in each local energy conservation domain are allocated and resided in the corresponding local memory physical address segment.

3. The method for improving the accuracy of multiphysics coupled simulation solutions according to claim 2, characterized in that, In S1, the step of generating the spatial distance attenuation weight matrix includes: Collect historical perturbation propagation datasets of the global coupling interface under different grid sizes, extract the pre-obtained optimal reference weight distribution as a supervision label, construct a loss function to iteratively update the spatial decay feature length constant, and train a mathematical model of boundary distance decay kernel function. Traverse any discrete grid node within the local energy conservation domain and calculate the shortest Euclidean distance scalar to the cross-domain boundary of the corresponding local energy conservation domain. Substitute this scalar into the boundary distance decay kernel function mathematical model to generate the corresponding dimensionless weight value. Construct the spatial distance decay weight matrix using all dimensionless weight values ​​as diagonal elements.

4. The method for improving the accuracy of multiphysics coupled simulation solutions according to claim 1, characterized in that, In S2, the step of extracting the numerical stiffness feature norm of the local energy conservation domain and the spatial gradient extrema of the feature physical quantities to calculate the tolerance threshold includes: Extract the system Jacobian matrix or stiffness matrix generated after discretizing the local energy conservation domain control equations, and perform Frobenius norm calculation to generate the numerical stiffness characteristic norm. In the multiphysics simulation variables, a preset master characteristic physical quantity is locked, the gradient vector of the physical quantity of the discrete grid node in three-dimensional space is calculated, and the gradient norm with the largest value inside the local energy conservation domain is selected as the spatial gradient extremum of the characteristic physical quantity. The numerical stiffness feature norm is multiplied by the spatial gradient extremum of the feature physical quantity and substituted into the empirical scaling constant obtained through data-driven training to unify the dimensions, thereby generating the tolerance threshold for local numerical stability.

5. The method for improving the accuracy of multiphysics coupled simulation solutions according to claim 1, characterized in that, In S3, the step of asynchronously advancing the source field solver and the target field solver and performing spatial interpolation mapping at the coupling communication point includes: Extract the minimum grid size and medium wave velocity of the source field solver and the target field solver respectively, and configure the source field time step and the target field time step to be an integer fraction of the preset macroscopic coupling communication period according to the preset Courant-Friedrich-Levy stability condition. Monitor the local advance clocks of the source field solver and the target field solver, and when both are extrapolated to the corresponding coupling communication point, trigger a synchronization barrier command to pause the time extrapolation process and enter the space data exchange stage. A spatial transformation matrix between the source field mesh and the target field mesh is constructed by calling the radial basis function interpolation algorithm. The spatial interpolation mapping result of nodal forces and displacement increments on the global coupling interface is calculated using the spatial transformation matrix.

6. The method for improving the accuracy of multiphysics coupled simulation solutions according to claim 1, characterized in that, In S4, the step of calculating the difference between the virtual work output by the source field and the virtual work received by the target field and accumulating them into a cumulative energy residual includes: Extract the source field nodal force vector output by the source field solver and the target field displacement increment vector output by the target field solver; The discrete virtual work transmitted from the source field physical domain to the global coupling interface is subtracted from the discrete virtual work received by the target field physical domain from the global coupling interface by the discrete inner product operation, and the discrete energy residual within the current preset macroscopic coupling communication cycle is calculated. The absolute value of the discrete energy residual is compared with the tolerance threshold. When the absolute value is less than the tolerance threshold, the discrete energy residual is algebraically summed with the historical residual value stored in the corresponding local energy residual accumulator, and the result is updated and written back to the local system to form the accumulated energy residual.

7. The method for improving the accuracy of multiphysics coupled simulation solutions according to claim 1, characterized in that, In S5, the step of determining the instability risk based on the accumulated energy residual and triggering an interrupt signal includes: In the established background monitoring thread, the absolute value of the accumulated energy residual is compared with the preset global tolerance limit in real time. When the absolute value exceeds the global tolerance limit, the running state transitions from the normal propulsion state to the alarm preparation state and the time window timer is started. The energy divergence rate is generated by calculating the average rate of change of the absolute value of the accumulated energy residual within a set sliding time window. The energy divergence rate is compared with the instability critical slope. When the energy divergence rate is greater than the instability critical slope, it is determined that the system has an instability risk, and the operating state transitions to the interrupt execution state, triggering the non-blocking type of interrupt signal.

8. The method for improving the accuracy of multiphysics coupled simulation solutions according to claim 7, characterized in that, The critical slope for instability is obtained through the following steps: The temporal energy fluctuation trajectories of different geometric configuration meshes running near the critical divergence boundary are extracted in advance, and the corresponding energy divergence rate features are extracted as training datasets. The state of whether or not a physical collapse is actually detected by the multiphysics simulation engine is defined as a system monitoring label; Based on the support vector machine classification algorithm, a classification hinge loss function is constructed between the predicted state and the true state. The decision boundary of the maximum margin hyperplane that distinguishes between stable and collapsed states is obtained by kernel function mapping and Lagrange multiplier optimization, which serves as the instability critical slope.

9. The method for improving the accuracy of multiphysics coupled simulation solutions according to claim 6, characterized in that, In S6, the step of converting the accumulated energy residual into an additional nodal force vector by combining the spatial distance attenuation weight matrix includes: Calculate the Euclidean norm of the displacement increment vector of each node in the target field and construct the node displacement weight coefficient. Multiply the coefficient with the corresponding dimensionless weight value in the spatial distance attenuation weight matrix to construct the composite spatial distribution weight. The composite spatial distribution weight is used for reverse mapping and multiplied with the displacement increment direction vector corresponding to each node to generate the basic direction force vector. The cumulative energy residual is divided by the sum of the virtual work contributions of all basic direction force vectors to obtain the energy equivalent scaling factor, which is multiplied into the basic direction force vector to generate the initial error compensation load. Extract the discrete stiffness matrix and mass matrix corresponding to the target field control equation to establish a generalized eigenvalue solution equation, iteratively obtain a low-order rigid body mode vector set to construct an orthogonal projection operator, and perform spatial filtering on the initial error compensation load through the orthogonal projection operator to remove invalid energy components and generate the additional nodal force vector.

10. The method for improving the accuracy of multiphysics coupled simulation solutions according to claim 9, characterized in that, In S6, the step of superimposing the external source terms into the target field solver includes: The generated additional node force vector is transformed into a dynamic source term function that decays according to an exponential function, wherein the time decay constant in the dynamic source term function is configured to be 0.5 to 2 times the preset macroscopic coupling communication period; The dynamic source term function is injected as an external load into the right-hand side of the non-homogeneous linear equation system of the target field solver to perform a source term reconstruction operation, and a synchronous release command is issued to unlock the synchronous barrier and restore the asynchronous time extrapolation process. After the action time window corresponding to the target field solver ends, the external load is removed to complete the energy smoothing correction, and the numerical data in the local energy residual accumulator is cleared to perform a cyclic iteration operation.

Images