A gpu-accelerated fluid-structure interaction simulation method and apparatus

Abstract

The application provides a GPU-accelerated fluid-structure interaction simulation method and device, and applies to the technical field of high-performance computing and industrial software crossing. The method comprises the following steps: in the matrix-free tree code radial basis function mapping, in each round of fluid-structure interaction iteration, the cross-GPU communication of external ghost data and the local near-far field calculation are deeply overlapped by using a multi-stream asynchronous strategy, and a spatially decoupled mixed precision scheme combined with Horner polynomial approximation is adopted to eliminate the calculation bottleneck of the logarithmic transcendental function ln in the thin plate spline kernel function; in this way, the GPU instruction pipeline stall is eliminated, the instruction-level floating point throughput is improved, the network communication delay under a super large-scale cluster is effectively hidden, and the "mapping wall" bottleneck is broken.

Description

Technical Field

[0001] This specification relates to one or more embodiments in the field of high-performance computing and industrial software intersection technology, and in particular to a GPU-accelerated fluid-structure interaction simulation method and apparatus. Background Technology

[0002] Fluid-structure interaction (FSI) simulation of nuclear reactors is central to nuclear reactor safety assessment and advanced design. In FSI simulations with ultra-large-scale grids (e.g., billions to tens of billions of grids), modern heterogeneous supercomputers face two major computational bottlenecks: one is the "mapping wall" problem, where data transfer between fluid-structure mismatched grid interfaces (such as global radial basis function (RBF) mapping) generates extremely high global communication overhead, and the large number of logarithmic transcendental functions (ln) contained in the underlying thin plate spline (TPS) kernel function can severely clog the GPU instruction pipeline.

[0003] Therefore, there is an urgent need for a GPU-accelerated fluid-structure interaction simulation method that can eliminate GPU instruction pipeline stalls, improve instruction-level floating-point throughput, effectively hide network communication latency in ultra-large-scale clusters, and break the "mapping wall" bottleneck. Summary of the Invention

[0004] This invention describes a GPU-accelerated fluid-structure interaction simulation method and apparatus, which can solve the above-mentioned technical problems.

[0005] According to the first aspect, a GPU-accelerated fluid-structure interaction simulation method is provided. The method includes: in each round of fluid-structure interaction iteration, performing interface mapping based on the original physical quantities of the interface output by the current solver in the previous iteration, and the first GPU instance performing mapping data preparation to obtain the payload; performing cross-node data exchange based on the payload through message passing communication to enable the current computing node to obtain the remote data required for interface mapping; recovering the remote data on the first GPU instance based on the exchanged payload to obtain the first calculation result corresponding to the external data; the current solver is a fluid solver or a structure solver; and during the message passing communication, the second GPU instance performs local near-field direct calculation and local far-field approximation calculation to obtain the second calculation result corresponding to the local data; accumulating the first calculation result and the second calculation result to obtain an output vector; solving a system of linear equations based on the first kernel function matrix and the output vector to obtain the radial basis function weight coefficient vector corresponding to the current iteration; the first kernel function matrix is ​​implemented in matrix-free form and obtained by substituting the distance between the corresponding source interface grid points into the thin plate spline kernel function; the thin plate spline kernel function is obtained by introducing static regularization parameters to reconstruct the kernel function and approximating it with Horner polynomials.

[0006] In some embodiments, the method further includes: in the multiphysics solution stage, constructing a bi-chain dependency graph for the heterogeneous resource requirements of the fluid solver and the underlying core operators of the structure solver; including: parsing the DGeMM task set of the fluid solver and the SpMV task set of the structure solver respectively, and constructing initial task dependency graphs respectively; grouping mutually independent tasks into parallel execution stages through topological sorting to form tensor core task chains and CUDA core task chains; each of the tensor core task chains and CUDA core task chains includes several execution stages, and the tasks in each execution stage are mutually independent; based on the cross tensor core task chains and CUDA core task chains, determining a cross-chain dependency edge set, the cross-chain dependency edge set including a first dependency pair and a second dependency pair, the first dependency pair including the last execution stage in the tensor core task chain and a unified synchronization point, the second dependency pair including the last execution stage in the CUDA core task chain and the unified synchronization point, the unified synchronization point being used to trigger synchronization to transition to the interface mapping stage.

[0007] In some embodiments, the method further includes: after constructing the initial task dependency graphs respectively, the method further includes: creating a Tensor Core Streaming Pool and a CUDA Core Streaming Pool; queuing independent tasks in each execution stage of the Tensor Core Task Chain to different streams in the Tensor Core Streaming Pool; queuing independent tasks in each execution stage of the CUDA Core Task Chain to different streams in the CUDA Core Streaming Pool; creating GPU events based on the dependency relationships of elements in the first dependency pair; recording the GPU events in the streams where the terminal tasks of the Tensor Core Task Chain and the CUDA Core Task Chain are located respectively; starting the next stage of the physics solving stage when the GPU events meet a preset state; simultaneously distributing thread blocks from the Tensor Core Streaming Pool and the CUDA Core Streaming Pool to the streaming multiprocessor (SM), and having the hardware instruction scheduler concurrently issue execution instructions to the Tensor Core and the CUDA Core according to the instruction mixing characteristics of different cores.

[0008] In some embodiments, the thin-plate spline kernel function is obtained by reconstructing the kernel function using a static regularization parameter and approximating it with a Horner polynomial; including:

[0009] Introducing static regularization parameters The kernel function of thin plate splines The equivalent reconstruction is:

[0010]

[0011] in To normalize the distance, it is pre-computed at compile time using constant folding. In near-field conditions Next, dynamic items Replace with a fifth-order Taylor expansion and rearrange into Horner polynomial form:

[0012]

[0013] In some embodiments, the second GPU instance performs local near-field direct computation and local far-field approximation computation, including: for the far-field approximation computation, using double-precision floating-point numbers (FP64) to prevent catastrophic cancellation during tree reduction; and for the near-field direct computation, using single-precision floating-point numbers (FP32) to double the theoretical floating-point throughput and reduce register pressure.

[0014] In some embodiments, the messaging communication includes:

[0015] Perform domain decomposition on the discrete mesh of fluid or structural domains with mesh mismatch;

[0016] Each GPU is responsible for a local region and constructs an octree structure for the grid points within that local region;

[0017] Each GPU obtains the spatial geometry of the remote local partition by exchanging the bounding boxes of the local octree root nodes, and generates a continuous one-dimensional communication array for recording message passing communication between nodes based on the following formula; the formula introduces a dynamic acceptance criterion of anisotropic penalty factor.

[0018]

[0019] in, The coordinates of the center of the local octree node. The center coordinates of the remote bounding box. Let be the radius of the bounding sphere of the local octree node. The spatial diagonal vector of the remote bounding box. Connect the vectors to the center ( The angle between the fluid-structure interaction interface and the principal deformation axis. An anisotropic penalty factor is introduced to address the elongated geometric features of the interface. This is the preset dynamic tolerance threshold.

[0020] In some embodiments, generating a continuous one-dimensional communication array for recording message passing communication between nodes includes: obtaining the spatial bounding boxes of all remote GPU local mesh partitions to determine the minimum spatial geometry of other partitions; each GPU is responsible for a local region and constructs an octree structure for the mesh points within that local region; for each obtained remote bounding box, the local GPU performs a depth-first search on its local octree; during the traversal, the tree nodes are dynamically classified using formula (3). When the currently traversed local tree node satisfies the above formula (3) If the local node is classified as far-field data, the aggregated unfolding moment of the local node is extracted, and the downward search of the branch is immediately terminated; if the leaf node still does not satisfy formula (3), it is classified as near-field data, and the precise spatial point coordinates contained in the leaf node are extracted. The extracted far-field node data and near-field point data are packaged and serialized into a one-dimensional communication array with contiguous memory; the one-dimensional communication array records the data load that the local GPU needs to send to each remote Rank, thereby transforming the irregular pointer-chasing tree structure into a regular memory layout adapted to asynchronous MPI transmission.

[0021] In some embodiments, a first GPU instance performs mapping data preparation to obtain a payload; based on the payload, performs cross-node data exchange via message passing communication to enable the current computing node to obtain the remote data required for interface mapping; based on the exchanged payload, the remote data is restored on the first GPU instance to obtain a first computation result corresponding to the external data; including: performing local moment calculation and communication data packaging, transmitting the payload to the host and immediately initiating non-blocking sparse point-to-point MPI communication; when the MPI communication is completed, the received external ghost data is returned, and the first GPU instance performs unpacking and ghost data evaluation; during the MPI communication, the second GPU instance sequentially performs local near-field interaction and local far-field evaluation, so that the computation of the second GPU instance and the MPI communication deeply overlap; and securely accumulates the first computation result of the first GPU instance and the second computation result of the second GPU instance into the global output vector through a write-after-write lock.

[0022] In some embodiments, the method further includes: updating the interface boundary conditions corresponding to the fluid solver or solid solver based on the radial basis function weight coefficient vector corresponding to the current iteration; determining whether the current iteration satisfies the preset convergence condition; if not, continuing to execute the next fluid-structure interaction iteration; if so, outputting the fluid-structure interaction simulation results.

[0023] According to a second aspect, a GPU-accelerated fluid-structure interaction simulation apparatus is provided, the apparatus comprising:

[0024] The first execution module is used to perform interface mapping based on the original physical quantities of the interface output by the current solver in the previous iteration in each round of fluid-structure interaction iteration. The first GPU instance performs mapping data preparation to obtain the payload. Based on the payload, cross-node data exchange is performed through message passing communication so that the current computing node obtains the remote data required for interface mapping. Based on the exchanged payload, the remote data is restored on the first GPU instance to obtain the first calculation result corresponding to the external data. The current solver is a fluid solver or a structure solver.

[0025] The second execution module is used to, during the message passing communication, have the second GPU instance perform local near-field direct computation and local far-field approximate computation to obtain a second computation result corresponding to local data; and to accumulate the first computation result and the second computation result to obtain an output vector.

[0026] The solution module is used to solve a system of linear equations based on the first kernel function matrix and the output vector to obtain the radial basis function weight coefficient vector corresponding to the current iteration. The first kernel function matrix is ​​implemented in a matrix-free form and is obtained by substituting the distance between the corresponding source interface grid points into the thin plate spline kernel function. The thin plate spline kernel function is obtained by introducing static regularization parameters to reconstruct the kernel function and approximating it with Horner polynomials.

[0027] In the methods and apparatus provided in the embodiments of this specification, a spatially decoupled hybrid precision scheme combined with Horner polynomial approximation replaces the near-field logarithmic transcendental function with five register-level fused multiply-accumulate FMA instructions. Simultaneously, a multi-stream asynchronous strategy is used to hide cross-node MPI communication latency. This completely eliminates GPU instruction pipeline stalls, significantly improving instruction-level floating-point throughput; effectively hides network communication latency in ultra-large-scale clusters, breaks the "mapping wall" bottleneck, and achieves efficient parallel scaling at the scale of tens of thousands of GPUs. Attached Figure Description

[0028] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0029] Figure 1 An architecture diagram of a GPU-accelerated fluid-structure interaction simulation method provided in this application embodiment;

[0030] Figure 2 A schematic diagram of a GPU-accelerated fluid-structure interaction simulation method provided in an embodiment of this application;

[0031] Figure 3 A schematic diagram of a continuous one-dimensional communication array constructed based on octree traversal, provided for an embodiment of this application;

[0032] Figure 4 An execution timing diagram of a multi-stream asynchronous and communication overlap mechanism in an MFTC-RBF mapping step provided in this application embodiment;

[0033] Figure 5This application provides a GPU-accelerated fluid-structure interaction simulation device. Detailed Implementation

[0034] The solution provided in this specification will now be described with reference to the accompanying drawings.

[0035] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions in the embodiments of this application will be described below with reference to the accompanying drawings.

[0036] In the description of the embodiments of this application, the words "exemplary," "for example," or "for instance" are used to indicate examples, illustrations, or explanations. Any embodiment or design described as "exemplary," "for example," or "for instance" in the embodiments of this application should not be construed as being more preferred or advantageous than other embodiments or designs. Specifically, the use of the words "exemplary," "for example," or "for instance" is intended to present the relevant concepts in a specific manner.

[0037] In the description of the embodiments of this application, the term "and / or" is merely a description of the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent: A existing alone, B existing alone, and A and B existing simultaneously. Furthermore, unless otherwise stated, the term "multiple" means two or more.

[0038] The terms “including,” “comprising,” “having,” and variations thereof all mean “including but not limited to,” unless otherwise specifically emphasized.

[0039] To address the above issues, this application provides a GPU-accelerated fluid-structure interaction (FSI) simulation method. Two concurrent streams are instantiated on the GPU, enabling asynchronous overlapping communication between the first and second GPU instances. The first GPU instance performs tree data preparation (local moment calculation and communication data packaging), transmits the payload to the CPU via Device-to-Host (D2H), and the host immediately initiates non-blocking sparse point-to-point (P2P) MPI communication. During MPI communication, the second GPU instance sequentially performs local near-field precise interaction and local far-field approximation evaluation locally. Since these calculations do not rely on external ghost data, their computation process deeply overlaps with network I / O, effectively hiding communication latency. After MPI communication is complete, the received external ghost data is transmitted back to the GPU via Host-to-Device (H2D), where the first GPU instance unpacks it and performs ghost data evaluation. Then, a Write-After-Write (WAW) locking mechanism is used to safely accumulate the computation results of the first and second GPU instances into the global output vector, and stream synchronization is performed.

[0040] Figure 1This is an architecture diagram of a GPU-accelerated fluid-structure interaction simulation method provided in an embodiment of this application. Figure 1 As shown, the architecture of the fluid-structure interaction simulation method includes a preprocessing module, a mapping module, and a DAS solution module.

[0041] The preprocessing module is primarily used to perform region decomposition on discrete meshes of fluid or structural domains with mesh mismatches. Each GPU is responsible for a local region and constructs a local octree structure for the mesh points within that region. Each GPU obtains the spatial geometry of the remote local partition by exchanging the bounding boxes of the root nodes of the local octrees. Subsequently, a depth-first search is performed on the local octrees, based on the Dynamic Acceptance Criterion (MAC) formula, which incorporates an anisotropy penalty factor.

[0042]

[0043] in, The coordinates of the center of the local octree node. The center coordinates of the remote bounding box. Let be the radius of the bounding sphere of the local octree node. The spatial diagonal vector of the remote bounding box. Connect the vectors to the center ( The angle between the fluid-structure interaction interface and the principal deformation axis. An anisotropic penalty factor is introduced to address the elongated geometric features of the interface. This is the preset dynamic tolerance threshold.

[0044] The mapping step module is primarily used to decouple global dense matrix-vector multiplication into local near-field exact calculation, local far-field approximation calculation, and external ghost data calculation in the Krylov subspace iteration of solving interface physical quantities. The mapping step module includes a communication computation overlap module and a mixed-precision and polynomial approximation module.

[0045] The communication-computation overlap module is primarily used to construct concurrent far-field streams (first GPU instance) and near-field streams (second GPU instance) on the GPU. The far-field stream is responsible for packaging data and initiating non-blocking sparse point-to-point MPI communication with the host; during this period, the near-field stream sequentially performs local near-field and far-field computations. After the MPI communication is completed, the far-field stream unpacks and evaluates the received ghost data, and finally securely accumulates the results of the two streams using a write-after-write lock, achieving deep overlap between communication and computation.

[0046] The mixed-precision and polynomial approximation module is primarily designed for the far-field computation, which involves a large number of high-order moments. Double precision (FP64) is used to prevent catastrophic cancellation; for near-field computation, single precision (FP32) is employed to double the throughput and reduce register pressure. Furthermore, a static regularization parameter ϵ is introduced to optimize the thin-plate spline TPS kernel function. Reconstruction, in near-field conditions ( The following is a fifth-order Taylor expansion and rearrangement into Horner polynomial form: This operation completely eliminates the log-transcendental function ln, which occupies the GPU Special Function Unit (SFU), and replaces it with five register-level fused multiply-accumulate (FMA) instructions.

[0047] The DAS solution step module is primarily used to address the heterogeneous resource requirements of the underlying core operators of the fluid solver and structure solver by constructing a bi-chain dependency graph based on fine-grained task dependency analysis. The DAS solution step module includes a task mapping and pipeline construction module and a lock-free synchronization and concurrent execution module.

[0048] The task mapping and pipeline construction module is mainly used to form tensor core task chains and CUDA core task chains through topological sorting, and to scan and generate cross-chain dependency edge sets. It queues computationally constrained large-batch dense matrix multiplication (DGeMM) tasks in the fluid solver to the tensor cores pool; and queues memory-access-constrained large-batch sparse matrix-vector multiplication (SpMV) tasks in the structure solver to the CUDA cores pool.

[0049] The lock-free synchronization and concurrent execution module is primarily used for cross-chain dependencies, leveraging lightweight GPU events to achieve non-blocking device-side synchronization between producer and consumer streams. The scheduler simultaneously distributes thread blocks from two stream pools to the streaming multiprocessor (SM). The hardware instruction scheduler, based on the instruction mixing characteristics of different kernels, enables the concurrent execution of two types of orthogonal task requirements within the same SM, maximizing hardware utilization.

[0050] Finally, assess whether the current coupling time step has reached the preset maximum time step. If not, proceed to the next time step and trigger the next dynamic iteration; otherwise, terminate the FSI simulation calculation.

[0051] Figure 2 This is a schematic diagram of a GPU-accelerated fluid-structure interaction simulation method provided in an embodiment of this application. Figure 2 As shown, the method includes at least the following steps:

[0052] S210. In each round of fluid-structure interaction iteration, interface mapping is performed based on the original physical quantities of the interface output by the current solver in the previous round. The first GPU instance performs mapping data preparation to obtain the payload. Based on the payload, cross-node data exchange is performed through message passing communication so that the current computing node obtains the remote data required for interface mapping. Based on the exchanged payload, the remote data is restored on the first GPU instance to obtain the first calculation result corresponding to the external data. The current solver is a fluid solver or a structure solver.

[0053] In one possible implementation, the current-side solver is the solver that is currently providing the interface data in this round of fluid-solid coupling iteration. If the current operation is a structure-to-fluid interface mapping, then the current-side solver is the solid solver.

[0054] In one possible implementation, such as Figure 4 As shown, the first GPU instance executes far-field flow. Part of the operation involves the second GPU instance performing near-field flow. Partial operations; far-field flow Some operations correspond to the acquisition and calculation of external data, near-field flow. Some operations correspond to calculations on local data. For example... Figure 2 As shown, one Rank is bound to one GPU; for example, Rank0 is bound to GPU0. Each GPU includes far-field flow. (First GPU instance) and near-field flow (Second GPU instance).

[0055] For example, during the execution of interface mapping computation, the first GPU instance needs to package some interface-related information and send it to other ranks, while also receiving similar information from other ranks. This packaged data is called the payload. For example, the payload is the set of relevant data that needs to be exchanged for performing interface mapping computation. Remote data comes from the payload sent by remote nodes.

[0056] In one possible implementation, prior to performing interface mapping, processing steps are executed, including communication topology and data structure construction. Communication topology and data structure construction includes:

[0057] S201. Perform region decomposition on the discrete mesh of the fluid or structural domains where the mesh does not match. Each GPU in the cluster first constructs a local octet structure based on the assigned local mesh points.

[0058] Figure 3 This is a schematic diagram illustrating the structure of a continuous one-dimensional communication array constructed based on octree traversal, as provided in an embodiment of this application. Figure 3 As shown, for example, each GPU in the cluster (e.g. Figure 3 The Rank0, Rank1, and Rank2 octrees in the diagram construct a local octree structure based on the assigned local grid points. After performing a depth-first search on the octree corresponding to Rank0, the black dots on the left side of the diagram represent the source points in the Rank0 local tree. After two-dimensional projection, nodes n1 (near) and n4 (near) that have near-field interactions with the target object (e.g., Rank1), and nodes n2 (far) and n3 (far) that have far-field interactions with the target object (e.g., Rank1) can be identified. Based on this near-field and far-field partitioning, Rank0 extracts and packages the near-field point data corresponding to the near-field nodes and the node aggregation data corresponding to the far-field nodes, and sends the data to be exchanged to Rank1, enabling Rank1 to perform near-field direct calculations and far-field approximate calculations for its target points.

[0059] S202. Each GPU obtains the spatial bounding box of the local mesh partition through a single NCCL All-gather set communication, thereby determining the minimum spatial geometric range of other partitions.

[0060] For example, the spatial bounding box is (Bounding Box = ).

[0061] S203. For each acquired remote bounding box, the local GPU performs a depth-first search (DFS) on its local octree. During the traversal, the tree nodes are dynamically classified using a dynamic acceptance criterion (MAC) that incorporates an anisotropic penalty factor.

[0062]

[0063] in, The coordinates of the center of the local octree node. The center coordinates of the remote bounding box. Let be the radius of the bounding sphere of the local octree node. The spatial diagonal vector of the remote bounding box. Connect the vectors to the center ( The angle between the fluid-structure interaction interface and the principal deformation axis. An anisotropic penalty factor is introduced to address the elongated geometric features of the interface. This is the preset dynamic tolerance threshold.

[0064] S204. If the currently traversed local tree node satisfies the above MAC formula (i.e.) If the data is not satisfied with the MAC formula when traversing to a leaf node, it is classified as "far-field data", the aggregation unfolding moment of the node is extracted, and the downward search of the branch is terminated immediately; if the MAC formula is not satisfied when traversing to a leaf node, it is classified as "near-field data", and the coordinates of the precise spatial points contained in the leaf node are extracted.

[0065] S205. Pack all extracted far-field node data and near-field point data into a contiguous one-dimensional communication array. This array explicitly records the data payload that the local GPU needs to send to each remote Rank, thereby transforming the irregular pointer-chasing tree structure into a regular memory layout adapted to asynchronous MPI transmission.

[0066] S220, and during the message passing communication, the second GPU instance performs local near-field direct calculation and local far-field approximate calculation to obtain a second calculation result corresponding to the local data; the first calculation result and the second calculation result are accumulated to obtain an output vector.

[0067] In one possible implementation, the core task in the FGMRES (Generalized Minimum Residual Method) Krylov subspace iteration for solving interface physical quantities is to compute the dense matrix-vector product q=Ap. This application's embodiment decouples and accelerates this process through the following steps:

[0068] In one possible implementation, at the start of each iteration, the local GPU needs to recalculate the moments of the entire local octree because the test vector p is updated.

[0069] First, the initial moments are calculated using a Taylor expansion of the leaf nodes with the latest weights. Then, multipole-to-multipole operators dependent on the geometric center offset Δr between nodes are used. Recursively accumulate child node moments from bottom to top Get parent node moments :

[0070]

[0071] In one possible implementation, two concurrent streams are instantiated on the GPU: the far-field stream and the far-field stream. and near-field flow . Figure 4 This is an execution timing diagram of a multi-stream asynchronous and communication overlap mechanism in an MFTC-RBF mapping step provided in an embodiment of this application. (See diagram for example.) Figure 4 As shown,

[0072] 1. First, tree data preparation (local moment calculation and communication data packaging) is performed. The payload is then transferred to the CPU via Device to Host (D2H), and the host immediately initiates non-blocking sparse point-to-point (P2P) message passing interface (MPI) communication.

[0073] 2. During MPI communication, Local near-field precise interaction and local far-field approximate evaluation are performed sequentially locally. Because these calculations do not rely on external ghost data, their computation process deeply overlaps with network I / O, effectively hiding communication latency.

[0074] 3. After MPI communication is complete, the received external ghost data is transmitted back to the GPU via Host to Device (H2D) communication. Unpack and perform ghost data evaluation.

[0075] 4. Finally, a write-after-write (WAW) locking mechanism is adopted to... and The calculation results are safely accumulated into the global output vector q, and stream synchronization is performed.

[0076] In one possible implementation, for far-field approximation calculations, due to the hierarchical aggregation of a large number of higher-order moments in the octree, double-precision floating-point numbers (FP64) are forced to prevent catastrophic cancellation during the tree reduction process; for near-field accurate calculations, since the range is strictly limited to local leaf nodes and their adjacent bounding boxes, single-precision floating-point numbers (FP32) are used for calculations, which doubles the theoretical floating-point throughput and reduces register pressure.

[0077] S230. Based on the first kernel function matrix and the output vector, solve the linear equation system to obtain the radial basis function weight coefficient vector corresponding to the current iteration; the first kernel function matrix is ​​implemented in matrix-free form and is obtained by substituting the distance between the corresponding source interface grid points into the thin plate spline kernel function; the thin plate spline kernel function is obtained by introducing static regularization parameters to reconstruct the kernel function and approximating it with Horner polynomials.

[0078] In one possible implementation, the TPS kernel function is used for thin plate splines. To address the problem of the log-transcendental function ln severely blocking GPU special function units (SFUs), a static regularization parameter is introduced. The kernel function is equivalently reconstructed as follows:

[0079]

[0080] in This is the normalized distance. It is pre-computed at compile time using constant folding. In near-field conditions Next, dynamic items Replace with a fifth-order Taylor expansion and rearrange into Horner polynomial form:

[0081]

[0082] This refactoring replaces the originally intensive ln calculations with five register-level fused multiply-accumulate FMA instructions, completely eliminating instruction pipeline stalls.

[0083] Traditional methods directly compute the global RBF matrix or use standard treecode, which suffers from communication congestion and performance bottlenecks due to the GPU's Special Function Unit (SFU) calculating the logarithmic function (ln). This application's embodiments introduce an anisotropic dynamic acceptance criterion to optimize the treecode topology and, through a spatially decoupled hybrid precision scheme combined with Horner polynomial approximation, replace the near-field logarithmic transcendental function with five register-level fused multiply-accumulate FMA instructions. Simultaneously, a multi-stream asynchronous strategy is used to hide cross-node MPI communication latency. This completely eliminates GPU instruction pipeline stalls, significantly improving instruction-level floating-point throughput; effectively hides network communication latency in ultra-large-scale clusters, breaks the "mapping wall" bottleneck, and achieves efficient parallel scaling at tens of thousands of GPUs.

[0084] In one possible implementation, to address the heterogeneous resource requirements of the fluid solver (computationally limited, with its core being large-scale DGeMM dense matrix multiplication) and the structure solver (memory-access limited, with its core being large-scale SpMV sparse matrix-vector multiplication), the following scheduling is performed:

[0085] S501 (Dual-chain dependency graph construction): Analyze the DGeMM task set of the fluid solver separately. SpMV task set for structure solvers Construct the initial task dependency graph and .

[0086] Specifically, the solution process is not a single giant matrix operation, but consists of hundreds or thousands of subtasks. Each DGeMM or SpMV operation is treated as a node, and the data dependencies between tasks (i.e., the order of execution) are treated as directed edges, thus constructing two directed acyclic graphs (DAGs). Among them, the fluid task (DGeMM) is computationally intensive and mapped to the GPU's Tensor Cores (TC); the structure task (SpMV) is memory intensive and mapped to regular CUDA cores (CUDA Cores, CC).

[0087] S502 (Topology Sorting and Cross-Chain Analysis): Groups independent tasks into parallel execution phases through topology sorting, forming a tensor core task chain. and CUDA core task chain Then, the cross-chain dependency edge set is defined. It is used to capture the dependencies between the end tasks of two solution chains and the unified synchronization point.

[0088] In one possible implementation, topological sorting extracts independent tasks in the graph that have no sequential dependencies and packages them into the same execution phase. For example, It includes the first set of fluid tasks that can be fully parallelized. By sorting, the originally chaotic tasks are organized into two independent pipelines that can be deployed in batches.

[0089] In one possible implementation, the aforementioned "two solution chains" refer to the fluid task chains mapped to the tensor core. With the structured task chain mapped to the CUDA core "End task" refers to the last execution phase in each of the two DAGs that has no successor node (i.e., the fluid chain). and structural chain The "unified synchronization point" represents the complete completion of the physical field calculations within the current iteration step; the "unified synchronization point" refers to the logical barrier used to ensure that the data on both sides of the fluid-structure interaction are ready before entering the next physical stage (i.e., the mesh mapping stage).

[0090] In one possible implementation, cross-chain dependency resolution and The structure is determined by the Jacobi coupling scheme; the fluid and structure solvers execute completely independently during the iteration phase, ensuring that intermediate dependencies across solvers are strictly empty. In terms of data structure, It is a set of directed edges, internally storing "producer-consumer" task tuples. Since there is no blocking in between, There are only two dependency pairs internally: and .in, and As producers, they generate the final flow field pressure and structural displacement, while the mesh mapping task... As consumers, they are located at a unified synchronization point. This allows two heterogeneous pipelines to proceed concurrently without any intermediate blockages, only passing through at the end. Trigger synchronization to safely transition to subsequent stages.

[0091] S503 (Pipeline Construction): Creates Tensor Core Stream Pool and CUDA core stream pool .Will Independent tasks in the queue to In different flows, Independent tasks in the queue to In different flows.

[0092] In one possible implementation, a "stream" in the GPU represents a queue of instructions executed sequentially, with different streams able to execute concurrently. The system ensures that fluid tasks and structural tasks are physically isolated from each other and do not interfere with each other by pre-allocating two dedicated stream pools.

[0093] S504 (Lockless Cross-Stream Synchronization): For... Each dependency in Create lightweight GPU events. In producer tasks... The event is recorded in the stream and in the consumer task. The device waits for the event in its stream, thus achieving non-blocking device-side synchronization and avoiding global barriers on the CPU side.

[0094] In one possible implementation, traditional methods typically require the CPU to issue a global synchronization instruction when handling multiphysics stage transitions (causing the entire GPU to idle and wait). This application uses the Jacobi format, eliminating intermediate dependencies in the solution stage, and this step does not require any intermediate synchronization intervention during the iteration process.

[0095] The embodiments of this application can achieve lock-free concurrency throughout the entire solution iteration. and All task flows run completely independently and at full capacity, without any intermediate pauses or blockages. Only when both pipelines have reached their end and triggered GPU events will the hardware automatically allow and safely transition to the mesh mapping stage, achieving true concurrency without intermediate blocking.

[0096] S505 (Concurrent Execution): The scheduler uses the LAUNCHALL instruction to send data from... and The thread blocks are simultaneously distributed to the GPU's streaming multiprocessor (SM). Due to the significantly different instruction mixing characteristics of the two types of cores, the hardware warp scheduler can concurrently issue execution instructions to Tensor Cores and CUDA Cores, maximizing resource utilization within a single SM.

[0097] In one possible implementation, thanks to the complete independence and lock-free event synchronization provided by the Jacobi format, computationally intensive fluid tasks and memory-intensive structural tasks break down the "serial alternation execution" barrier. Both types of tasks run simultaneously at the physical level on the same SM (Single Memory Module) of the same GPU, without any intermediate blocking throughout the iteration phase. Tensor cores handle intensive multiply-accumulate operations, while CUDA cores and memory bandwidth handle sparse memory access operations, maximizing resource utilization by squeezing the throughput of heterogeneous hardware to its limit.

[0098] In one possible implementation, the next iteration begins after the current iteration is completed. This includes:

[0099] S601: At the current time step After the solution and interface mapping are completed, the time controller evaluates whether the preset maximum simulation time step has been reached.

[0100] S602: If not achieved, proceed to the next time step. Update the fluid-structure interface boundary conditions using the mapping results of the current step, and trigger the next dynamic iteration (return to step S220).

[0101] S603: If the maximum time step has been reached, output the final flow field and structural deformation simulation results and terminate the entire FSI simulation calculation.

[0102] Traditional coupling schemes execute fluid and structure solvers serially or in isolation on different computing nodes, resulting in alternating idle hardware resources. This invention constructs a bi-chain dependency graph based on fine-grained task dependency analysis, directing computationally intensive tasks of the fluid solver to Tensor Cores and memory-intensive tasks of the structure solver to CUDA Cores, and utilizing lightweight GPU events to achieve lock-free cross-stream synchronization. This overcomes the isolation limitations at the solver level, enabling concurrent execution of two types of tasks with orthogonal resource requirements within the same streaming multiprocessor (SM). It transforms the previously accumulated solution time into overlapping execution, maximizing the utilization of heterogeneous hardware resources within the GPU and significantly shortening the overall solution time for fluid-structure interaction (FSI) problems.

[0103] Figure 5 This application provides a GPU-accelerated fluid-structure interaction simulation device. For example... Figure 5 As shown, the device includes:

[0104] The first execution module is used to perform interface mapping based on the original physical quantities of the interface output by the current solver in the previous iteration in each round of fluid-structure interaction iteration. The first GPU instance performs mapping data preparation to obtain the payload. Based on the payload, cross-node data exchange is performed through message passing communication so that the current computing node obtains the remote data required for interface mapping. Based on the exchanged payload, the remote data is restored on the first GPU instance to obtain the first calculation result corresponding to the external data. The current solver is a fluid solver or a structure solver.

[0105] The second execution module is used to, during the message passing communication, have the second GPU instance perform local near-field direct computation and local far-field approximate computation to obtain a second computation result corresponding to local data; and to accumulate the first computation result and the second computation result to obtain an output vector.

[0106] The solution module is used to solve a system of linear equations based on the first kernel function matrix and the output vector to obtain the radial basis function weight coefficient vector corresponding to the current iteration. The first kernel function matrix is ​​implemented in a matrix-free form and is obtained by substituting the distance between the corresponding source interface grid points into the thin plate spline kernel function. The thin plate spline kernel function is obtained by introducing static regularization parameters to reconstruct the kernel function and approximating it with Horner polynomials.

[0107] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made on the basis of the technical solution of the present invention should be included within the scope of protection of the present invention.

Claims

1. A GPU-accelerated fluid-structure interaction simulation method, characterized in that, The method includes: In each iteration of fluid-structure interaction, interface mapping is performed based on the original physical quantities of the interface output by the current solver in the previous iteration. The first GPU instance prepares the mapping data to obtain the payload. Based on the payload, cross-node data exchange is performed through message passing communication so that the current computing node obtains the remote data required for interface mapping. Based on the exchanged payload, the remote data is recovered on the first GPU instance to obtain the first calculation result corresponding to the external data. The current solver is a fluid solver or a structure solver. During the message passing communication, the second GPU instance performs local near-field direct computation and local far-field approximate computation to obtain a second computation result corresponding to the local data; the first computation result and the second computation result are accumulated to obtain an output vector; Based on the first kernel function matrix and the output vector, a system of linear equations is solved to obtain the radial basis function weight coefficient vector corresponding to the current iteration. The first kernel function matrix is ​​implemented in matrix-free form and is obtained by substituting the distance between the corresponding source interface grid points into the thin plate spline kernel function. The thin plate spline kernel function is obtained by reconstructing the kernel function and approximating it with Horner polynomials by introducing static regularization parameters. The thin plate spline kernel function is obtained by reconstructing the kernel function and approximating it with Horner polynomials by introducing static regularization parameters, including: introducing static regularization parameters. The kernel function of thin plate splines The equivalent reconstruction is: in To normalize the distance, it is pre-computed at compile time using constant folding. In near-field conditions Next, dynamic items Replace with a fifth-order Taylor expansion and rearrange into Horner polynomial form: 。 2. The method according to claim 1, characterized in that, The method further includes: In the multiphysics solution phase, a bi-chain dependency graph is constructed to address the heterogeneous resource requirements of the underlying core operators of the fluid solver and the structure solver; including: The DGeMM task set at the bottom layer of the fluid solver and the SpMV task set at the bottom layer of the structure solver are analyzed separately, and the initial task dependency graphs are constructed respectively. Independent tasks are grouped into parallel execution stages through topological sorting, forming the Tensor Core Task Chain and the CUDA Core Task Chain; each of the Tensor Core Task Chain and the CUDA Core Task Chain includes several execution stages, and the tasks in each execution stage are independent of each other. Based on the Tensor Core Task Chain and the CUDA Core Task Chain, a cross-chain dependency edge set is determined. The cross-chain dependency edge set includes a first dependency pair and a second dependency pair. The first dependency pair includes the last execution stage in the Tensor Core Task Chain and a unified synchronization point. The second dependency pair includes the last execution stage in the CUDA Core Task Chain and the unified synchronization point. The unified synchronization point is used to trigger synchronization to transition to the interface mapping stage.

3. The method according to claim 2, characterized in that, The method further includes: after constructing the initial task dependency graphs respectively, the method further includes: Create the Tensor Core Streaming Pool and the CUDA Core Streaming Pool; The independent tasks in each execution stage of the tensor core task chain are queued into different streams of the tensor core stream pool; The independent tasks in each execution stage of the CUDA core task chain are queued into different streams in the CUDA core stream pool; Based on the dependency relationships of the elements in the first dependency pair, create a GPU event; The GPU events are recorded in the streams containing the terminal tasks of the Tensor Core Task Chain and the CUDA Core Task Chain, respectively; when the GPU events meet the preset state, the next stage of the physics solving stage is started. Thread blocks from the Tensor Core Stream Pool and the CUDA Core Stream Pool are simultaneously distributed to the Streaming Multiprocessor (SM). The hardware instruction scheduler then concurrently issues execution instructions to the Tensor Core and the CUDA Core based on the instruction mixing characteristics of different cores.

4. The method according to claim 1, characterized in that, The second GPU instance performs local near-field direct computation and local far-field approximation computation, including: For the far-field approximation calculation, double-precision floating-point numbers are used to prevent catastrophic cancellation during the tree reduction process; For the near-field direct calculation, single-precision floating-point numbers are used for calculation, which doubles the theoretical floating-point throughput and reduces register pressure.

5. The method according to any one of claims 1 to 4, characterized in that, The message passing communication includes: Perform domain decomposition on the discrete mesh of fluid or structural domains with mesh mismatch; Each GPU is responsible for a local region and constructs an octree structure for the grid points within that local region; Each GPU obtains the spatial geometry of the remote local partition by exchanging the bounding boxes of the local octree root nodes, and generates a continuous one-dimensional communication array for recording message passing communication between nodes based on the following formula; the formula introduces a dynamic acceptance criterion of anisotropic penalty factor. in, The coordinates of the center of the local octree node. The center coordinates of the remote bounding box. Let be the radius of the bounding sphere of the local octree node. The spatial diagonal vector of the remote bounding box. Connect the vectors to the center ( The angle between the fluid-structure interaction interface and the principal deformation axis. An anisotropic penalty factor is introduced to address the elongated geometric features of the interface. This is the preset dynamic tolerance threshold.

6. The method according to claim 5, characterized in that, The generation of a continuous one-dimensional communication array for recording message passing communication between nodes includes: Obtain the spatial bounding boxes of all remote GPU local mesh partitions to determine the minimum spatial geometry of other partitions; Each GPU is responsible for a local region and constructs an octree structure for the grid points within that local region; For each acquired remote bounding box, the local GPU performs a depth-first search on its local octree; During the traversal, formula (3) is used to dynamically classify the tree nodes. When the currently traversed local tree node satisfies the above formula (3) If the local tree node is classified as far-field data, the aggregation unfolding moment of the local tree node is extracted, and the downward search of the branch is terminated immediately; if the leaf node is still not satisfied with formula (3), it is classified as near-field data, and the precise spatial coordinates contained in the leaf node are extracted. The extracted far-field node data and near-field point data are packaged and serialized into a one-dimensional communication array with contiguous memory. The one-dimensional communication array records the data payload that the local GPU needs to send to each remote Rank, thereby transforming the irregular pointer-chasing tree structure into a regular memory layout adapted to asynchronous MPI transmission.

7. The method according to claim 1, characterized in that, The first GPU instance performs mapping data preparation to obtain a payload; based on the payload, it performs cross-node data exchange through message passing communication so that the current computing node obtains the remote data required for interface mapping. Based on the swapped payload, remote data is recovered on the first GPU instance to obtain the first calculation result corresponding to the external data; include: Perform local moment calculation and communication data packaging, transmit the payload to the host and immediately initiate non-blocking sparse point-to-point MPI communication; after the MPI communication is completed, the received external ghost data is transmitted back and unpacked and ghost data evaluated by the first GPU instance; During the MPI communication, the second GPU instance sequentially performs local near-field interaction and local far-field evaluation, so that the computation of the second GPU instance is deeply overlapped with the MPI communication. The first computation result of the first GPU instance and the second computation result of the second GPU instance are safely accumulated into the global output vector by using a write-after-write lock.

8. The method according to claim 1, characterized in that, The method further includes: Based on the radial basis function weight coefficient vector corresponding to the current iteration, update the interface boundary conditions corresponding to the fluid solver or solid solver. Determine whether the current iteration meets the preset convergence condition; if not, continue to execute the next round of fluid-structure interaction iteration; if so, output the fluid-structure interaction simulation results.

9. A GPU-accelerated fluid-structure interaction simulation device, characterized in that, The device includes: The first execution module is used to perform interface mapping based on the original physical quantities of the interface output by the current solver in the previous iteration in each round of fluid-structure interaction iteration. The first GPU instance performs mapping data preparation to obtain the payload. Based on the payload, cross-node data exchange is performed through message passing communication so that the current computing node obtains the remote data required for interface mapping. Based on the exchanged payload, the remote data is restored on the first GPU instance to obtain the first calculation result corresponding to the external data. The current solver is a fluid solver or a structure solver. The second execution module is used to, during the message passing communication, have the second GPU instance perform local near-field direct computation and local far-field approximate computation to obtain a second computation result corresponding to local data; and to accumulate the first computation result and the second computation result to obtain an output vector. The solution module is used to solve a system of linear equations based on the first kernel function matrix and the output vector to obtain the radial basis function weight coefficient vector corresponding to the current iteration. The first kernel function matrix is ​​implemented in matrix-free form and is obtained by substituting the distance between the corresponding source interface grid points into the thin plate spline kernel function. The thin plate spline kernel function is obtained by introducing static regularization parameters to reconstruct the kernel function and approximating it with Horner polynomials. The thin plate spline kernel function is obtained by introducing static regularization parameters to reconstruct the kernel function and approximating it with Horner polynomials; including: introducing static regularization parameters. The kernel function of thin plate splines The equivalent reconstruction is: in To normalize the distance, it is pre-computed at compile time using constant folding. In near-field conditions Next, dynamic items Replace with a fifth-order Taylor expansion and rearrange into Horner polynomial form: 。

Images