Multi-GPU cooperative parallel computing method for fast environmental wind field simulation

Abstract

The present application relates to the technical field of wind field numerical simulation and graphics processor (GPU) parallel computing, and discloses a multi-GPU cooperative parallel computing method for fast environmental wind field simulation, which comprises the following steps in sequence: GPU resource configuration, three-dimensional computing domain X direction strip decomposition, persistent video memory allocation, initial data distribution, divergence parallel computing, multi-GPU cooperative successive over-relaxation iteration solving, velocity field correction and result assembly. The method adopts global coordinate red-black coloring to guarantee the coloring consistency across subdomains, combines asynchronous ghost layer data exchange and hierarchical two-stage error reduction, supports peer-to-peer direct access and host relay dual-path communication, and eliminates repeated allocation overhead through the persistent video memory mechanism. The present application can realize multi-GPU efficient cooperative solving of Poisson equation, significantly improves the computing efficiency and resource utilization rate of large-scale urban wind field simulation, and is suitable for various hardware topological structures.

Description

Technical Field

[0001] This invention relates to the field of wind field numerical simulation and parallel computing technology of graphics processing units (GPUs), specifically to a multi-GPU collaborative parallel computing method for rapid environmental wind field simulation. Background Technology

[0002] Rapid Environmental Simulation (QES) wind field solutions are based on a variational analysis framework. By solving the Poisson equation with Lagrange multipliers, urban wind fields that satisfy mass conservation constraints are obtained, which has significant engineering value in fields such as urban planning, atmospheric environmental assessment, and building wind load analysis. As the resolution of urban wind field simulations continues to increase and the computational domain scale continues to expand, the memory capacity and computing power of a single GPU can no longer support the needs of large-scale, high-precision simulations, resulting in a significant increase in solution time, which has become the core bottleneck limiting the efficiency of wind field simulations.

[0003] Existing multi-GPU parallel wind field solution schemes mostly employ a heterogeneous architecture combining message passing interfaces and GPUs. Data interaction between devices relies on host memory relay, resulting in high communication latency and large transmission overhead. Most schemes use dynamic memory allocation, leading to frequent memory allocation and deallocation when repeatedly calling the solver. Furthermore, existing red-black sorting iterations generally use subdomain local coordinate coloring rules, which easily result in adjacent grid points of the same color across subdomain boundaries, disrupting the dependencies in iterative calculations and compromising the accuracy of the solution results.

[0004] Furthermore, existing multi-GPU parallel methods lack hardware topology adaptation capabilities and cannot automatically switch communication paths based on device peer-to-peer access capabilities, resulting in insufficient compatibility under different hardware configurations. Global error convergence judgment often adopts cross-device atomic reduction methods, which have high synchronization overhead between devices, further limiting the parallel acceleration efficiency and making it difficult to meet the application requirements of QES wind field simulation for efficient, universal, and highly reliable parallel solutions. Summary of the Invention

[0005] To address the aforementioned problems, this invention provides a multi-GPU collaborative parallel computing method for fast environmental wind field simulation, comprising the following steps:

[0006] Step S1, System Initialization and GPU Resource Configuration: Complete GPU hardware resource discovery, topology detection and peer access configuration, build peer access reachability matrix, enable bidirectional direct access for devices that support peer direct memory access, and automatically fall back to the host relay path if enabling fails. Step S2, 3D computational domain spatial decomposition: Divide the global 3D structured computational domain into subdomains matching the number of GPUs along the X direction, maintain the global complete scale in the Y and Z directions, add ghost layers with a width of 1 to the non-global boundary side of the subdomains, and configure global-local coordinate mapping and data extraction and assembly rules. Step S3, Persistent Memory Resource Allocation: Allocate independent memory resources and computation flow to each GPU. A persistent allocation mechanism is adopted. Memory allocation is performed when the solver function is called for the first time. Subsequent calls directly reuse the allocated memory space. All memory resources are released uniformly when the solver is destructed. Step S4, Initial Data Parallel Distribution: Extract subdomain data containing the ghost layer from the global wind field data and asynchronously upload it to the corresponding GPU's video memory to complete the initial data initialization; Step S5: Parallel calculation of initial wind field divergence: Each GPU calculates the divergence array of the right-hand side of the Poisson equation in parallel based on the initial velocity field of its local subdomain. Step S6: Multi-GPU collaborative successive over-relaxation iterative solution: A global coordinate red-black coloring mechanism is adopted to ensure coloring consistency across subdomains. Red grid points and black grid points are alternately updated half-step by half-step. After each half-step iteration, asynchronous ghost layer data exchange is performed to synchronize boundary data. Hierarchical global error reduction is used for convergence judgment. Iteration continues until the global error is less than the convergence tolerance or the maximum number of iterations is reached. Step S7, Parallel correction of velocity field: Based on the converged Lagrange multiplier field, the initial velocity field is corrected in parallel to obtain the final wind field that satisfies the divergence-free condition; Step S8, Result Collection and Global Assembly: Copy back the subdomain calculation results of each GPU, extract the data of the non-overlapping regions of the subdomain, and assemble them to obtain the complete global wind field.

[0007] Preferably, in step S1, the total number of GPUs available in the system is enumerated by the device quantity acquisition function built into the unified computing device architecture, all GPU devices are traversed to query the peer-to-peer direct memory access capability, and the query results are stored as a Boolean two-dimensional peer-to-peer access reachability matrix; during solver destructing, the peer-to-peer access disable function is called to release resources for all device pairs that have enabled peer-to-peer access.

[0008] Preferably, in step S2, when dividing in the X direction, the number of basic units and the number of remaining units of a single GPU are calculated, and one more unit is allocated to each of the subdomains corresponding to the remaining number to achieve balanced distribution of computing load; the width of the ghost layer matches the first-order neighborhood data requirements of the 7-point template of the successive over-relaxation solver, and is used to store the boundary data of adjacent subdomains, so that each GPU can complete the iterative calculation of the boundary grid points locally.

[0009] Preferably, in step S3, the memory resources of each GPU are encapsulated as a single-device data structure, which includes a body-centered array, a face-centered array, an auxiliary array, and an independent unified computing device architecture stream handle; wherein the body-centered array includes a successive over-relaxation boundary condition coefficient array, a divergence array, a Lagrange multiplier array, an array of the previous iteration values, and a grid cell flag array, and the face-centered array includes velocity component arrays in the x, y, and z directions.

[0010] Preferably, in step S6, the formula for determining the global coordinate red-black coloring is as follows:

[0011] in For local x-direction index, This represents the global x-coordinate value corresponding to the left boundary of the subdomain. Grid points are defined as follows: grid points that meet this condition are red grid points, and grid points with a color offset parameter of 1 correspond to black grid points, thus avoiding the problem of adjacent grid points with the same color due to local coordinate coloring of subdomains.

[0012] Preferably, in step S6, during asynchronous ghost layer data exchange, all data in column X of the subdomain boundary are first... Elements are packaged into a continuous one-dimensional send buffer; devices that support peer-to-peer direct memory access transmit data through direct asynchronous copying between video memory, and ensure the timing correctness of data transmission and unpacking through cross-stream event synchronization; devices that do not support peer-to-peer access automatically switch to a host relay path.

[0013] Preferably, in step S6, the hierarchical global error reduction adopts a two-level reduction architecture: each GPU calculates the maximum absolute error in its local subdomain through atomic comparison and exchange operations to complete device-level local error reduction; after synchronizing the computing flow of all GPUs, the local errors of each device are copied back to the host, and the maximum value is taken as the global maximum error to complete host-level global error aggregation.

[0014] Preferably, in step S6, the relaxation factor for each successive over-relaxation iteration is 1.78; after each complete iteration, each GPU performs the boundary condition application operation in parallel, applying Neumann boundary conditions at the bottom boundary of the computation domain.

[0015] Preferably, in step S7, the initial velocity components in the x, y, and z directions are corrected based on the Euler-Lagrange equation, and the velocity values ​​in the solid / terrain unit are set to zero to obtain the final wind field that satisfies the divergence-free condition.

[0016] Preferably, in step S4, the data upload operations of each subdomain are executed in parallel through an independent unified computing device architecture stream, synchronizing the computing streams of all devices to ensure that all data transmissions are completed before proceeding to subsequent computing steps.

[0017] Compared with the prior art, the beneficial effects of the present invention are as follows: One-dimensional striped domain decomposition in the X direction is used to achieve balanced distribution of computational load. A global coordinate red-black coloring mechanism ensures the correctness of cross-subdomain iterative dependencies. Asynchronous ghost layer data exchange reduces communication waiting overhead. Hierarchical two-level error reduction avoids cross-device atomic operation losses. The computational scale can be linearly expanded with the number of GPUs, which greatly shortens the solution time of large-scale high-resolution wind fields.

[0018] It adopts a persistent memory allocation mechanism, and after the initial memory allocation is completed, subsequent calls can directly reuse the memory, eliminating the performance overhead of repeated allocation and release; it supports GPU peer-to-peer direct memory access to reduce communication latency, and automatically falls back to the host relay path when access fails, which can adapt to different hardware topologies and ensure general availability in multi-device scenarios.

[0019] The ghost layer width matches the neighborhood data requirements of the 7-point iterative template, ensuring accurate iterative calculation of boundary grid points; the velocity field correction is completed based on the variational analysis framework, strictly satisfying the divergence-free mass conservation constraint; the overall process is adapted to the QES wind field solver architecture and can be directly integrated into engineering scenarios for rapid environmental simulation of urban wind fields. Attached Figure Description

[0020] Figure 1 This is an overall flowchart of the method of the present invention.

[0021] Figure 2 This is a schematic diagram of the two-dimensional grid red-black sorting iteration of the present invention. Detailed Implementation

[0022] This invention discloses a multi-GPU (Graphics Processing Unit) collaborative parallel computing method for wind field simulation in QES (Fast Environmental Simulation). Applied to the wind field solver of a fast environmental simulation system, it solves urban wind fields satisfying mass conservation constraints based on a variational analysis framework. The core method employs a successive over-relaxation red-black sorting iterative method to solve the Poisson equation with Lagrange multipliers. This method achieves efficient collaborative parallel computing across multiple GPUs by decomposing a three-dimensional structured computational domain into one-dimensional stripes along the X-direction, combined with persistent memory management, global coordinate red-black coloring, asynchronous ghost layer data exchange, and hierarchical global error reduction techniques.

[0023] The overall solution process of this method includes eight stages: system initialization and graphics processor resource configuration, 3D computational domain spatial decomposition, persistent memory allocation, initial data distribution, divergence parallel computation, multi-graphics processor collaborative successive over-relaxation iterative solution, velocity field parallel correction, result collection, and global assembly. (Refer to...) Figure 1 The specific implementation methods for each step are as follows: Step 1: System Initialization and Graphics Processor Resource Configuration This step completes the discovery, topology detection, and peer-to-peer access configuration of graphics processor hardware resources, providing the hardware foundation for multi-graphics processor collaboration. The specific implementation process is as follows: 1. Solver Instantiation: The user specifies the solver type as a multi-GPU solver (type code 5) via the command line and can choose to specify the number of GPUs to use; the solver factory creation method instantiates a multi-GPU solver object based on the type code. If the user specifies the number of GPUs, the smaller value between the user-specified value and the actual number of GPUs available in the system is used; if the user-input parameter is 1, all available GPUs are used automatically.

[0024] 2. Device Enumeration: The total number of usable graphics processors in the system is obtained through the built-in device count function of the unified computing device architecture, thus completing device resource discovery.

[0025] 3. Construction of Peer Reachability Matrix: Iterate through all graphics processor device pairs, call the device peer access query function to query whether each pair of graphics processors supports peer direct memory access, and store the query results as a Boolean two-dimensional peer access reachability matrix; a true element in the matrix indicates that the corresponding device pair supports peer direct memory access, and a false element indicates that it does not.

[0026] 4. Peer Access Enable and Rollback: Peer access is enabled for all device pairs whose values ​​in the reach matrix are true. The device peer access enable function is called to activate bidirectional direct memory access. If enabling fails (e.g., due to system security policy restrictions), the corresponding matrix position is reset to false, and subsequent communication is automatically switched to the host relay path, ensuring availability under any device topology.

[0027] 5. Resource release mechanism: During solver destructing, the device peer access disable function is called for all device pairs with peer access enabled to release peer access resources.

[0028] Step 2: 3D computational domain spatial decomposition and subdomain configuration: This step divides the global 3D structured computation domain into subdomains matching the number of graphics processors, and configures the boundary ghost layer and coordinate mapping rules. The specific implementation process is as follows: 1. One-dimensional strip partitioning in the X direction: The global computational domain is uniformly divided along the X direction into... Subdomains ( (This refers to the number of graphics processors participating in the computation), each subdomain maintains a globally complete scale in the Y and Z directions. Let the global domain have a total of [number of processors] in the X direction. Each face node corresponds to Individual grid cells refer to the number of basic cells allocated to a single graphics processor. for:

[0029] Number of units remaining after allocation for:

[0030] The remaining number of subdomains are each allocated the number of basic units plus 1 unit, and the remaining subdomains are each allocated the number of basic units, thus achieving a balanced distribution of computational load.

[0031] 2. Ghost Layer Configuration: For the non-global boundary side of each subdomain, a ghost layer with a width of 1 is added. The specific rule is: if a neighboring subdomain exists to the left of the subdomain, add one ghost cell to the left end of the local X direction; if a neighboring subdomain exists to the right of the subdomain, add one ghost cell to the right end of the local X direction. The width of the ghost layer matches the first-order neighborhood data requirements of the successive over-relaxation solver's 7-point template (center point and 6 neighboring points ±x, ±y, ±z), and is used to store the boundary data of adjacent subdomains, enabling each graphics processor to complete the iterative calculation of boundary grid points locally.

[0032] 3. Subdomain Information Definition: The configuration information of each subdomain is described by the subdomain information structure, which includes: the assigned graphics processor device number, the subdomain's X / Y offset in the global domain, the local staggered mesh dimension including the ghost layer, the ghost layer width, the left and right neighbor subdomain indices (a value of 1 indicates that the corresponding side is the global physical boundary), and the local index calculation auxiliary method, including the local cell index calculation method and the local surface index calculation method, used for the mutual conversion between global coordinates and local coordinates.

[0033] 4. Data Extraction and Assembly Rules: Define subdomain cell data extraction methods and subdomain surface data extraction methods to extract subdomain data containing ghost layer overlapping regions from the global data array; define a global surface data assembly method to write the calculation results of each subdomain back to the global array. This method ensures that only non-overlapping self-owned data is written back by calculating the start and end indices of the subdomain's own region, thus avoiding duplicate overwriting.

[0034] Step 3: Persistent allocation of video memory resources: This step allocates independent video memory resources and computational flow to each graphics processor, and uses a persistent allocation mechanism to eliminate the overhead of duplicate allocation. The specific implementation process is as follows: 1. Persistent allocation control: Set a persistent allocation flag to control that video memory is allocated only when the solver function is called for the first time. Subsequent calls directly reuse the allocated video memory space. All video memory resources are released uniformly when the solver is destructed.

[0035] 2. Single-device memory architecture configuration: Define a single-device data structure for each graphics processor, encapsulating all the data and resources required for computation in a single subdomain, specifically including: Body center array: The array of successively relaxed boundary condition coefficients, including the arrays of body center coefficients e, f, g, h, m, and n, the divergence array, the array of Lagrange multipliers and the array of the previous iteration value, and the array of grid cell flags; Face-centered arrays: x-direction velocity component array, y-direction velocity component array, z-direction velocity component array; Auxiliary array: Vertical grid spacing array, local error scalar; Control resources: Independent unified computing device architecture stream handles.

[0036] 3. Memory and Stream Creation: Traverse each subdomain, switch to the corresponding graphics processor device in turn (call the function to set the current device), create an independent unified computing device architecture stream (call the stream creation function), allocate memory for all the above volume-centered and face-centered arrays (call the memory allocation function), and set the persistent allocation flag to true after completion.

[0037] Step 4: Initial Data Parallel Distribution and Initialization: This step distributes the global initial data to the subdomain video memory of each graphics processor, completing the data preparation before computation. The specific implementation process is as follows: 1. Subdomain Data Extraction: For each subdomain, the extraction method of the domain decomposer is called to extract the subdomain data block containing the ghost layer from the global wind field general data structure, including the weight coefficients of adjacent grid points in the positive x-direction. Weighting coefficients of adjacent grid points in the negative x-direction Weighting coefficients of adjacent grid points in the positive y-direction Weight coefficients of adjacent grid points in the negative y-direction Weighting coefficients of adjacent grid points in the positive z-direction The weight coefficient of adjacent grid points in the negative z direction The array includes a grid cell label array, a vertical grid spacing array, and initial velocity components at the face center in the horizontal, vertical, and longitudinal directions.

[0038] 2. Asynchronous data upload: The asynchronous memory copy function is called to asynchronously upload the data of each subdomain to the memory of the corresponding graphics processor. The upload operation of each device is executed in parallel through the independent unified computing device architecture. At the same time, the asynchronous memory clearing function is used to clear and initialize the Lagrange multiplier array with the Lagrange multiplier array of the previous step.

[0039] 3. Stream synchronization: Synchronize the unified computing device architecture stream of all graphics processors to ensure that all data transmissions are completed before proceeding to subsequent computing steps.

[0040] Step 5: Parallel calculation of initial wind field divergence: This step involves parallel calculation of the divergence R of the right-hand side of the Poisson equation based on the initial velocity field, preparing for subsequent iterative solutions.

[0041] The rapid environmental simulation wind field solution is based on a variational analysis framework, minimizing the objective functional. for:

[0042] In the formula: These represent the velocity components of the final wind field in the x, y, and z directions, respectively. These represent the velocity components of the initial wind field in the x, y, and z directions, respectively. These are the weighting coefficients for the x, y, and z directions, respectively. For Lagrange multipliers; Let be the divergence of the velocity field; This represents the total volume of the three-dimensional computational domain.

[0043] The variational derivation of the above functional yields the Poisson equation satisfied by the Lagrange multipliers:

[0044] , and These are the horizontal, vertical, and axial coordinates, respectively. In the formula The divergence of the initial wind field is calculated in the following discrete form:

[0045] In the formula: These are the grid cell indices for the x, y, and z directions, respectively; These are the grid step sizes in the x, y, and z directions, respectively; and Grid points The initial x-velocity at the interface between the positive and negative x-directions; and Grid points The initial y-velocity at the interface between the positive and negative y-directions; and Grid points The initial z-velocity at the interface between the positive and negative z directions.

[0046] The specific implementation of this step is as follows: the divergence calculation kernel function is launched in parallel on each graphics processor, and each device executes in parallel through an independent unified computing device architecture to calculate the divergence array of the right-hand side of the Poisson equation based on the initial velocity field of the local subdomain.

[0047] Step 6: Multi-GPU collaborative successive over-relaxation iterative solution: This step is the core of the method. It solves the Poisson equation using a red-black sorting successive over-relaxation iterative method, combined with global coordinate coloring, asynchronous ghost layer swapping, and layered error reduction to achieve multi-device collaborative iterative computation. The iterative loop is repeated until the global error is less than the convergence tolerance or the maximum number of iterations is reached. (Refer to...) Figure 2 The left image corresponds to the red grid point iteration step, where all red grid points marked 1 have no computational dependency and can be updated synchronously in parallel. The right image corresponds to the black grid point iteration step, where all black grid points marked 2 can be updated synchronously in parallel. These two half-steps are executed alternately to complete a full round of successive over-relaxation iterations, which is the core parallelization foundation for parallel solving of the Poisson equation. The specific sub-steps of this process are as follows: 6.1 Save the old values ​​from the iteration: Each graphics processor executes the Lagrange multiplier storage kernel function in parallel, copying the values ​​of the current Lagrange multiplier array to the previous Lagrange multiplier array for subsequent iteration error calculation.

[0048] 6.2 Parallel Update of Red Half-Step Grid Points: This step completes the Lagrange multiplier update for all red grid points, and uses a global coordinate red-black coloring mechanism to ensure coloring consistency across subdomains.

[0049] The discrete calculation formula for successive over-relaxation iterations is as follows:

[0050] In the formula: For grid points The Lagrange multiplier iteration value at the location; The successive over-relaxation relaxation factor is set to 1.78. These are two computational domain constants, where Anisotropic scaling corresponding to the y-direction, Anisotropic scaling corresponding to the z-direction; For grid points The divergence value at; Adjacent points in the positive x-direction Lagrange multipliers; Adjacent points in the negative x-direction Lagrange multipliers; Adjacent points in the positive y-direction Lagrange multipliers; Adjacent points in the negative y-direction Lagrange multipliers; Adjacent points in the positive z-direction Lagrange multipliers; Adjacent points in the negative z direction Lagrange multipliers.

[0051] Specific implementation process: Each graphics processor executes the successive over-relaxed red-black sorting kernel function in parallel (with the shading offset parameter set to 0). Within the kernel function, each thread recovers its local 3D coordinates based on a one-dimensional thread index, calculates the global X coordinates using the subdomain global X offset, and determines red-black shading using the following formula:

[0052] For local x-direction index; This represents the global x-coordinate value corresponding to the left boundary of the subdomain; it indicates the equality judgment.

[0053] Grid points that meet this condition are designated as red grid points. After skipping the ghost layer region, the thread executes the above successive over-relaxation iterative formula on the red grid points to complete the red half-step Lagrange multiplier update.

[0054] This global coordinate coloring mechanism avoids the problem of adjacent boundaries with the same color caused by local coordinate coloring in subdomains, ensures the consistency of chessboard coloring across subdomains, and ensures the correctness of computation in successive over-relaxation iterations.

[0055] 6.3 First asynchronous data exchange at the ghost layer: This step exchanges the updated red grid boundary data to the ghost layer of the neighboring subdomains, providing neighborhood data for the black half-step update. The specific implementation process is as follows: 1. Data Packaging: Call the ghost-level data packaging kernel function to pack all data in column X of the subdomain boundary. Elements are packed into a contiguous one-dimensional send buffer. Each thread in the unified computing device architecture processes one. The source data address of an element is calculated as follows:

[0056] Offset to the starting column. This represents the number of units in the local X direction. The number of units in the local Y direction is the total number of units in the subdomain X direction, and the number of units in the local Y direction is the total number of units in the subdomain Y direction.

[0057] 2. Peer-to-peer asynchronous transmission: For device pairs that support peer-to-peer direct memory access, the peer-to-peer asynchronous video memory copy function is called to achieve direct asynchronous data copying between video memory. This transmission does not pass through the host memory and achieves low-latency interaction through high-speed interconnect links or fast peripheral component interconnect buses.

[0058] 3. Cross-stream event synchronization: On the source device's unified computing device architecture stream, call the event recording function to record the event and mark the peer-to-peer transmission as complete; on the target device's unified computing device architecture stream, call the stream wait event function to wait for the event, ensuring that the unpacking operation is only started after the data transmission is complete, avoiding data races.

[0059] 4. Data unpacking: Call the ghost-level data unpacking kernel function to distribute the data from the continuous receive buffer back to the corresponding ghost column positions in the Lagrange multiplier array of the target device, thus completing the boundary data update.

[0060] 5. Non-peer fallback path: When the device does not support peer-to-peer direct memory access, it automatically switches to the host relay path: first, the data is copied from the source device to the host memory through the device-to-host mode, and then the data is copied from the host memory to the target device through the host-to-device mode, ensuring full topology compatibility.

[0061] 6.4 Parallel Update of Black Half-Step Grid Points: Each graphics processor executes the successive over-relaxed red-black sorting kernel function in parallel (with the shading offset parameter set to 1). The execution flow is symmetrical to the red half-step, completing the Lagrange multiplier update for all black grid points.

[0062] 6.5 Second Asynchronous Data Exchange in the Ghost Layer: Repeat the above ghost layer exchange process to transmit the updated black grid boundary data to the corresponding ghost layer region of the neighboring subdomain, completing one full iteration of boundary data synchronization.

[0063] 6.6 Parallel application of boundary conditions: Each graphics processor executes the von Neumann boundary condition application kernel function in parallel, applying the von Neumann boundary condition at the bottom boundary of the computation domain (k=0).

[0064] 6.7 Hierarchical Error Reduction and Convergence Judgment: A two-level reduction architecture is adopted to achieve global error aggregation across devices. The specific implementation process is as follows: 1. Device-level local error reduction: Each graphics processor resets its local local error scalar to 0, starts the error calculation kernel function, calculates the maximum absolute error in the local subdomain through atomic comparison and exchange operations, and completes the error reduction within a single device.

[0065] 2. Host-level global error aggregation: After synchronizing the unified computing device architecture flow of all graphics processors, the local error values ​​of each device are copied back to the host, and the maximum value of all local errors is taken as the global maximum error on the host.

[0066] 3. Convergence judgment: Compare the global maximum error with the convergence tolerance. If it is less than the convergence tolerance, exit the iteration loop; otherwise, continue to the next iteration.

[0067] This hierarchical reduction design avoids performance loss from cross-device atomic operations while ensuring the global correctness of convergence judgments.

[0068] Step 7, Parallel correction of velocity field: After iterative convergence, each graphics processor executes the final velocity field calculation kernel function in parallel. Based on the converged Lagrange multiplier field, the initial velocity field is corrected using the Euler-Lagrange equation to obtain the final wind field that satisfies the divergence-free condition. The correction formula is as follows:

[0069]

[0070]

[0071] The symbols in the formula have the same meaning as described above. Simultaneously, the kernel function sets the velocity values ​​in the solid / terrain cells to zero, completing the velocity field correction.

[0072] Step 8: Results Collection and Global Data Assembly This step aggregates the subdomain calculation results of each graphics processor into global wind field data. The specific implementation process is as follows: 1. Subdomain result copy: For each subdomain, the velocity component arrays in the x, y, and z directions in the device's video memory are copied back to the host temporary buffer using the video memory copy function.

[0073] 2. Global Data Assembly: Call the global surface data assembly method to extract the velocity data of each subdomain's own region (excluding ghost layer overlap areas), write it into a global array, assemble it to obtain a complete global velocity field, and return it to the caller.

Claims

1. A multi-GPU collaborative parallel computing method for rapid environmental wind field simulation, characterized in that, Includes the following steps: Step S1, System Initialization and GPU Resource Configuration: Complete GPU hardware resource discovery, topology detection and peer access configuration, build peer access reachability matrix, enable bidirectional direct access for devices that support peer direct memory access, and automatically fall back to the host relay path if enabling fails. Step S2, 3D computational domain spatial decomposition: Divide the global 3D structured computational domain into subdomains matching the number of GPUs along the X direction, maintain the global complete scale in the Y and Z directions, add ghost layers with a width of 1 to the non-global boundary side of the subdomains, and configure global-local coordinate mapping and data extraction and assembly rules. Step S3, Persistent Memory Resource Allocation: Allocate independent memory resources and computation flow to each GPU. A persistent allocation mechanism is adopted. Memory allocation is performed when the solver function is called for the first time. Subsequent calls directly reuse the allocated memory space. All memory resources are released uniformly when the solver is destructed. Step S4, Initial Data Parallel Distribution: Extract subdomain data containing the ghost layer from the global wind field data and asynchronously upload it to the corresponding GPU's video memory to complete the initial data initialization; Step S5: Parallel calculation of initial wind field divergence: Each GPU calculates the divergence array of the right-hand side of the Poisson equation in parallel based on the initial velocity field of its local subdomain. Step S6: Multi-GPU collaborative successive over-relaxation iterative solution: A global coordinate red-black coloring mechanism is adopted to ensure coloring consistency across subdomains. Red grid points and black grid points are alternately updated half-step by half-step. After each half-step iteration, asynchronous ghost layer data exchange is performed to synchronize boundary data. Hierarchical global error reduction is used for convergence judgment. Iteration continues until the global error is less than the convergence tolerance or the maximum number of iterations is reached. Step S7, Parallel correction of velocity field: Based on the converged Lagrange multiplier field, the initial velocity field is corrected in parallel to obtain the final wind field that satisfies the divergence-free condition; Step S8, Result Collection and Global Assembly: Copy back the subdomain calculation results of each GPU, extract the data of the non-overlapping regions of the subdomain, and assemble them to obtain the complete global wind field; In step S6, the formula for determining the red-black coloring of global coordinates is as follows: in For local x-direction index, This represents the global x-coordinate value corresponding to the left boundary of the subdomain. Grid points are defined as follows: grid points that meet this condition are red grid points, and grid points with a color offset parameter of 1 correspond to black grid points, thus avoiding the problem of adjacent grid points with the same color due to local coordinate coloring of subdomains.

2. The multi-GPU collaborative parallel computing method for rapid environmental wind field simulation according to claim 1, characterized in that, In step S1, the total number of GPUs available in the system is enumerated by the device quantity acquisition function built into the unified computing device architecture. All GPU devices are traversed to query the peer-to-peer direct memory access capability, and the query results are stored as a Boolean two-dimensional peer-to-peer access reachability matrix. When the solver is destructed, the peer-to-peer access disable function is called to release resources for all device pairs that have enabled peer-to-peer access.

3. The multi-GPU collaborative parallel computing method for rapid environmental wind field simulation according to claim 1, characterized in that, In step S2, when dividing in the X direction, the number of basic units and the number of remaining units of a single GPU are calculated. One more unit is allocated to each of the subdomains corresponding to the remaining number of units to achieve balanced distribution of computing load. The width of the ghost layer matches the first-order neighborhood data requirements of the 7-point template of the successive over-relaxation solver and is used to store the boundary data of adjacent subdomains, so that each GPU can complete the iterative calculation of the boundary grid points locally.

4. The multi-GPU collaborative parallel computing method for rapid environmental wind field simulation according to claim 1, characterized in that, In step S3, the memory resources of each GPU are encapsulated into a single-device data structure, which includes a body-centered array, a face-centered array, an auxiliary array, and an independent unified computing device architecture stream handle. The body-centered array includes a successive over-relaxation boundary condition coefficient array, a divergence array, a Lagrange multiplier array, an array of the previous iteration values, and a grid cell flag array. The face-centered array includes velocity component arrays in the x, y, and z directions.

5. The multi-GPU collaborative parallel computing method for fast environmental wind field simulation according to claim 1, characterized in that, In step S6, during asynchronous ghost layer data exchange, all data in column X of the subdomain boundary is first... Elements are packed into a contiguous one-dimensional send buffer; Devices that support peer-to-peer direct memory access transfer data through asynchronous copying between video memory, ensuring the timing correctness of data transmission and unpacking through cross-stream event synchronization; devices that do not support peer-to-peer access automatically switch to a host relay path.

6. The multi-GPU collaborative parallel computing method for fast environmental wind field simulation according to claim 1, characterized in that, In step S6, the hierarchical global error reduction adopts a two-level reduction architecture: each GPU calculates the maximum absolute error in its local subdomain through atomic comparison and exchange operations to complete device-level local error reduction; after synchronizing the computing flow of all GPUs, the local errors of each device are copied back to the host, and the maximum value is taken as the global maximum error to complete host-level global error aggregation.

7. The multi-GPU collaborative parallel computing method for fast environmental wind field simulation according to claim 1, characterized in that, In step S6, the relaxation factor for each successive over-relaxation iteration is 1.78; after each complete iteration, each GPU performs the boundary condition application operation in parallel, applying Neumann boundary conditions at the bottom boundary of the computation domain.

8. The multi-GPU collaborative parallel computing method for rapid environmental wind field simulation according to claim 1, characterized in that, In step S7, the initial velocity components in the x, y, and z directions are corrected based on the Euler-Lagrange equations, and the velocity values ​​in the solid / terrain units are set to zero to obtain the final wind field that satisfies the divergence-free condition.

9. The multi-GPU collaborative parallel computing method for fast environmental wind field simulation according to claim 1, characterized in that, In step S4, the data upload operations of each subdomain are executed in parallel through an independent unified computing device architecture stream. The computing streams of all devices are synchronized to ensure that all data transmissions are completed before proceeding to the subsequent computing steps.

Images