An openmp-based cpu / multi-gpu heterogeneous parallel acceleration three-dimensional cantilever beam topology optimization method and system

Abstract

The application belongs to the technical field of topology optimization, and discloses a CPU / multi-GPU heterogeneous parallel acceleration three-dimensional cantilever beam topology optimization method and system based on OpenMP. The method comprises the following steps: a subdivided grid design domain of a three-dimensional cantilever beam to be optimized is divided into multiple subdomains, and each subdomain is subjected to topology optimization; the optimized subdomains are spliced to obtain an optimized three-dimensional cantilever beam design domain; each subdomain is respectively optimized; S1 i=1, the residual of the current subdomain is calculated; S2 the residual of the current subdomain is updated, and the current subdomain is updated accordingly; S3 i=1+1, when i is less than or equal to a preset threshold N, step S2 is returned; otherwise, the residual of the Nth layer subdomain is solved; S4 each layer subdomain is reprocessed, the residual of the first layer is calculated, when the residual does not meet the preset threshold, step S2 is returned, otherwise, S5 is entered; S5 the material density of each grid unit is updated to update the first layer subdomain. Through the application, the parallel efficiency of topology optimization solving is improved.

Description

Technical Field

[0001] This invention belongs to the field of topology optimization technology, and more specifically, relates to a method and system for topology optimization of three-dimensional cantilever beams based on OpenMP and CPU / multi-GPU heterogeneous parallel acceleration. Background Technology

[0002] Topology optimization is a structural design method that seeks the optimal material distribution within a given design region based on specific physical constraints (such as volume and stress). It is widely used in lightweight design for high-end equipment in aerospace, automotive manufacturing, and other industries. With the increasing demands for design precision in industry, the mesh size of three-dimensional finite element analysis for typical load-bearing components such as three-dimensional cantilever beams is growing exponentially, often reaching tens or even hundreds of millions of degrees of freedom. Traditional CPU-based serial or MPI parallel solvers face unacceptably long computation times and extremely high distributed cluster hardware costs when handling such ultra-large-scale problems.

[0003] In recent years, graphics processing units (GPUs), with their massive concurrent execution threads, powerful floating-point arithmetic capabilities, and extremely high global memory bandwidth, have become the mainstream hardware platform for accelerating the solution of ultra-large-scale 3D cantilever beam topology optimization problems. However, for typical desktop single-node multi-GPU computing environments, existing parallel topology optimization schemes for 3D cantilever beams still have significant limitations in terms of hardware-software collaboration, facing the following four major technical bottlenecks: (1) "Memory Wall" Constraints and Bottlenecks in Ultra-Large-Scale Solving: In traditional topology optimization processes, the explicit assembly and storage of stiffness matrices is the main component of memory consumption. As the mesh size increases, the number of non-zero elements in the global sparse matrix increases dramatically. When a single consumer-grade or professional-grade GPU (e.g., the 24GB physical memory limit) processes models with hundreds of millions of degrees of freedom (e.g., more than 340 million degrees of freedom), it is very easy for the computation to crash directly due to memory overflow (OOM). Current general-purpose iterative solvers lack adaptive strategies for hybrid operators for different mesh resolutions, making it difficult to completely eliminate matrix storage in fine mesh layers with ultra-large degrees of freedom. This has become a physical computation ceiling that single-node hardware cannot overcome.

[0004] (2) Inefficient data layout and wasted memory bandwidth: The finite element template (Stencil) operator in topology optimization has low arithmetic strength and is a typical bandwidth-constrained application. Its performance limit is highly dependent on memory access efficiency. Most existing finite element parallel code uses the traditional array of structures (AoS) format, which physically stores the displacement degrees of freedom of a single node in the X, Y, and Z directions. However, in the single instruction multithreaded (SIMT) architecture of GPUs, this stepping access leads to severe non-merged memory access, generating a large number of redundant memory transactions, resulting in a large amount of extremely expensive high-speed memory bandwidth being idled and wasted. In addition, if the mesh size of the master order dimension does not meet the transaction byte alignment constraints of the underlying hardware, it will also cause cache line invalidation and doubled transmission latency.

[0005] (3) Lack of cross-card communication latency masking mechanism: To overcome the limitations of single-card memory, multi-GPU collaborative computing is an inevitable trend. However, when multi-grid algorithms (such as MGCG) move between different resolution levels, they need to perform constraints, interpolation, and halo layer boundary data exchange extremely frequently. In environments lacking direct hardware bridging and relying solely on the traditional PCIe bus, tiny data packet transmissions will face extremely high bus startup latency. Especially in small to medium-sized grids or extremely coarse grid layers, the extremely short internal core computation time cannot mask the high cross-card synchronization overhead, resulting in communication time accounting for more than 50%, which means that multi-GPU systems not only cannot achieve the ideal linear speedup, but may even experience a negative speedup phenomenon that is slower than single-card serial processing.

[0006] (4) The challenge of distributed consistency handling of physical boundaries: In large-scale parallel computing, in order to achieve load balancing, the global physical design domain is usually decomposed into multiple local subdomains along a certain coordinate axis (such as the X-axis). This geometric division directly leads to the fragmentation of the originally unified global physical boundary conditions (such as Neumann external load boundaries and Dirichlet fixed constraint boundaries). In matrix-free iterative solvers that abandon the global stiffness matrix, due to the lack of a unified coordinate system, there is still a lack of a standard and atomic efficient processing mechanism for how to apply these local load effects and constraint reactions distributed on different GPUs in parallel, consistently and without conflict on the device side with extremely low kernel startup overhead.

[0007] In summary, existing parallel topology optimization methods lack deep joint optimization from underlying memory layout and asynchronous pipelined communication to physical layer decomposition. There is an urgent need for a new CPU / multi-GPU heterogeneous parallel acceleration method that can completely break through the limitations of single-machine video memory, fully extract the peak computing power of heterogeneous hardware, and completely mask the latency of PCIe cross-card communication. Summary of the Invention

[0008] To address the aforementioned deficiencies or improvement needs of existing technologies, this invention provides a CPU / multi-GPU heterogeneous parallel acceleration method and system for three-dimensional cantilever beam topology optimization based on OpenMP, solving the problems of slow convergence and computational bottleneck in ultra-large-scale three-dimensional topology optimization.

[0009] To achieve the above objectives, according to one aspect of the present invention, a CPU / multi-GPU heterogeneous parallel acceleration method for topology optimization of three-dimensional cantilever beams based on OpenMP is provided, the method comprising the following steps: The subdivided mesh design domain of the 3D cantilever beam to be optimized is divided into multiple subdomains. Topology optimization is performed on each subdomain. The optimized subdomains are then joined together to obtain the optimized 3D cantilever beam design domain. The optimization of each subdomain is carried out according to the following steps: S1 i=1, take the grid cells of the current subdomain as the i-th layer, and calculate the residuals of each grid cell in the current subdomain; S2 updates the residuals of each grid cell in the current subdomain by multiplying the preset constraint operator with the residual, and updates each grid cell according to the updated residuals, thereby updating the current subdomain; S3 i=1+1, when i is less than or equal to the preset threshold N, return to step S2; otherwise, solve the residual of each grid cell in the Nth subdomain; S4 j=N, use the residuals of each grid cell in the j-th subdomain to reprocess the subdomain. This reprocessing includes interpolation, error correction and post-smoothing. j=j-1, until the reprocessing of all subdomains is completed. Calculate the residuals of the state equations of the grid cells in the first subdomain. If the residuals do not meet the preset threshold, return to step S2. Otherwise, proceed to S5. S5 Calculate the sensitivity of the material density of the first-layer subdomain mesh cell to compliance; update the material density of each mesh cell according to the volume constraint of the mesh cell and the sensitivity, and use the updated material density to update the first-layer subdomain mesh cell. The updated first-layer subdomain mesh cell is the optimized subdomain.

[0010] More preferably, after dividing the subdivided mesh design domain of the three-dimensional cantilever beam to be optimized into multiple subdomains, each subdomain is divided into an aura layer and a non-aura layer. Then, each subdomain is mapped to the GPU, and topology optimization is performed on each subdomain in the GPU. In step S3, when solving the state equation of the Nth subdomain, the design variables of the Nth layer need to be passed to the CPU connected to the GPU, and the residuals of each mesh cell of the Nth subdomain are solved in the CPU.

[0011] More preferably, in step S1, the residual of each grid cell is obtained by first solving the state equation of the subdomain to obtain the displacement of each grid cell, and then using the displacement to solve for the residual.

[0012] More preferably, the state equation of the subdomain is as follows: KU=F in, The subdomain stiffness matrix is ​​related to the material distribution. and These are the subdomain displacement vector and the subdomain external load vector, respectively.

[0013] More preferably, the formula for calculating the residual is as follows:

[0014] in, This is the residual.

[0015] More preferably, when calculating the state equation of a subdomain, for subdomains of layers 1 to N-1, the data of the aura layers in adjacent subdomains need to be exchanged in the CPU, and then the state equation of the current subdomain is calculated using the data of the aura layers in adjacent subdomains and the non-aura layers in the current subdomain; for subdomain of layer N, the data of the aura layers and non-aura layers in both adjacent subdomains and the current subdomain need to be transferred to the CPU for joint calculation, so as to obtain the state equation of the current subdomain.

[0016] More preferably, in step S4, the interpolation, error correction, and post-smoothing processes employ trilinear interpolation, direct summation, or damped Jacobian algorithms, respectively.

[0017] More preferably, in step S5, the formula for updating the material density of each grid cell based on the volume constraint of the grid cell and the sensitivity is as follows:

[0018] in, and These represent the relative material densities before and after the mesh cell update. The minimum density lower bound set for the singularity of the stiffness matrix. For the movement limit, The iteration factor is jointly determined by the objective function sensitivity and the volume constraint Lagrange multiplier. This is the numerical damping coefficient.

[0019] More preferably, in step S5, the formula for calculating the sensitivity of the material density of the first-layer subdomain mesh element to compliance is as follows:

[0020] in, The overall compliance of the structure, As a penalty factor, The material density of the current unit. This represents the nodal displacement vector of the current element. The material density is The element foundation stiffness matrix at that time.

[0021] According to another aspect of the present invention, a CPU / multi-GPU heterogeneous parallel acceleration system for three-dimensional cantilever beam topology optimization based on OpenMP is provided. The system includes an actuator for executing the aforementioned CPU / multi-GPU heterogeneous parallel acceleration method for three-dimensional cantilever beam topology optimization based on OpenMP.

[0022] In summary, the technical solutions conceived by this invention have the following beneficial effects compared with the prior art: 1. This invention divides the design domain to be optimized into multiple subdomains and updates each subdomain layer by layer using residuals, followed by a secondary update using material density. Existing technologies for ultra-large-scale 3D topology optimization (such as meshes with tens or hundreds of millions of nodes) require enormous amounts of memory / GPU memory to store the physical state field (displacement, residuals, etc.) and perform massive floating-point operations, which is insufficient for a single GPU's memory capacity. Therefore, the overall cantilever beam 3D mesh design domain must be divided into multiple subdomains along a certain dimension (along the X-axis), and each subdomain is loaded into different GPU memory locations, thus solving the slow convergence and computational bottleneck problems in ultra-large-scale 3D topology optimization.

[0023] 2. This invention maps each subdomain to multiple GPUs, thus freeing the problem size from the memory bottleneck of a single graphics card, and achieving a balanced distribution of memory pressure and computational load. Based on this architecture, the algorithm can be easily scaled to 4 cards, 8 cards, or even distributed clusters, thereby supporting ultra-high resolution (ultra-large scale) topology optimization structure design. By establishing a halo layer, adjacent GPUs can send and receive residual and displacement data of the boundary regions between each iteration (i.e., boundary data synchronization). This halo exchange strategy ensures that the results of multi-GPU parallel computation are completely equivalent to the results of computation using a single graphics card with ultra-large memory at the mathematical and physical mechanics level, without causing boundary stress concentration or loss of physical quantities due to region partitioning.

[0024] 3. This invention exchanges data from the aura layers in adjacent subdomains within the CPU and passes the design variables of the Nth layer to the CPU connected to that subdomain. The CPU then solves for the residuals of each mesh cell in the Nth subdomain. This design eliminates the need for top-tier supercomputing architectures in the parallel algorithm, significantly lowering the hardware threshold for ultra-large-scale topology optimization. This solves the computational bottlenecks in ultra-large-scale 3D topology optimization caused by the dramatic increase in degrees of freedom (tens of millions), such as "global matrix memory overflow," "limited multi-GPU computing power expansion," and "slow convergence of highly ill-conditioned systems."

[0025] 4. This invention transmits data to the CPU for joint computation at the coarsest layer (Nth layer) of the multi-grid, instead of continuing the iteration on the GPU. This avoids the negative speedup caused by frequent small data packet synchronization when multiple GPUs are cooperating, greatly masking communication latency and maximizing the heterogeneous parallel computing power of a single node with multiple GPUs.

[0026] 5. This invention employs The residual calculation mode utilizes constraint operators and interpolation operators at each mesh level to transmit error information, minimizing the memory consumption of a single topology optimization iteration. This enables ordinary node devices to overcome physical memory bottlenecks and smoothly run three-dimensional cantilever beam topology optimization with hundreds of millions or even larger degrees of freedom. It solves the problem that traditional rigid matrix explicit assembly consumes massive amounts of memory and is prone to overflow (OOM), thus avoiding the storage of huge global sparse matrices.

[0027] 6. This invention exchanges the aura layer data of adjacent subdomains through the CPU and performs joint state equation calculation with the non-aura layer data of the current subdomain. It uses physical grid partitioning to break the original global physical boundary conditions, avoiding the physical field after iteration from breaking at the partition surface. By abandoning the global unified coordinate system, it ensures the continuity and consistency of the physical field assembled after local update through precise overlap and synchronization of the aura layers, thereby ensuring the final accuracy of topology optimization. Attached Figure Description

[0028] Figure 1 This is a schematic diagram of one-dimensional region decomposition and aurora layer communication along the X-axis constructed according to a preferred embodiment of the present invention.

[0029] Figure 2 It is a multi-GPU collaborative multi-mesh inter-layer data transmission and physical boundary consistency mapping diagram constructed according to a preferred embodiment of the present invention.

[0030] Figure 3 This is a flowchart of the CPU / multi-GPU heterogeneous parallel multi-grid V-Cycle algorithm based on OpenMP, constructed according to a preferred embodiment of the present invention.

[0031] Figure 4 This is a schematic diagram of the topology optimization results of a three-dimensional cantilever beam under different mesh sizes using dual-GPU heterogeneous acceleration, constructed according to a preferred embodiment of the present invention.

[0032] Figure 5 This is the final configuration diagram of a high-resolution topology optimization of a three-dimensional cantilever beam with a scale of 512*256*256 under dual-GPU heterogeneous acceleration, constructed according to a preferred embodiment of the present invention. Detailed Implementation

[0033] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention. Furthermore, the technical features involved in the various embodiments of this invention described below can be combined with each other as long as they do not conflict with each other.

[0034] The global mesh design domain of the 3D cantilever beam to be optimized is divided and assigned to multiple GPU subdomains. In each topology optimization iteration, OpenMP is used to coordinate boundary communication between the subdomains to maintain the continuity of the global physics field, and each subdomain is optimized in parallel according to the following multi-mesh preprocessing conjugate gradient (MGCG) steps: S1: Within the outer conjugate gradient (CG) iterative framework, let The fine mesh of the current subdomain is taken as the first... Layer. The current layer is computed in parallel by the GPU. The state equation residuals of each grid cell in the layer are calculated, and pre-smoothing is performed on the residual field.

[0035] like Figure 1 As shown, after dividing the subdivided mesh design domain of the three-dimensional cantilever beam to be optimized into multiple subdomains, each subdomain is divided into an aura layer and a non-aura layer. Then, each subdomain is mapped to the GPU. Next, topology optimization is performed on each subdomain in the GPU. In step S3, when solving the state equation of the Nth subdomain, the design variables of the Nth layer need to be passed to the CPU connected to the GPU. The residuals of each mesh cell of the Nth subdomain are solved in the CPU.

[0036] This invention divides the subdivided mesh design domain of the 3D cantilever beam to be optimized into multiple subdomains, and further divides each subdomain into an aura layer and a non-aura layer, then maps each subdomain to a GPU. Since calculating the force and deformation of a mesh cell depends on the data of its surrounding neighboring cells, when the entire physical domain is divided and distributed to different GPUs, the meshes located at the edges of the subdomains will lose their original neighboring mesh data. Therefore, an overlapping data buffer layer, the "aura layer," needs to be artificially extended outward at the subdomain boundaries; while the internal regions that do not need to exchange data with other GPUs are called the "non-aura layer."

[0037] The residuals of each grid cell are obtained by first solving the state equations of the subdomain to obtain the displacements of each grid cell, and then using the displacements to solve for the residuals. The state equations of the subdomain are as follows: KU=F in, The subdomain stiffness matrix is ​​related to the material distribution. and These are the subdomain displacement vector and the subdomain external load vector, respectively.

[0038] The formula for calculating the residual is as follows:

[0039] in, This is the residual.

[0040] For any subfield Divide its node set into internal node sets Light Linkage Point Set Among them, optical nodes refer to points where the physical coordinates are located at... Within the range, but its state updates depend on neighboring subdomains. ( The nodes of the data. To ensure the numerical equivalence between parallel computing and serial algorithms, the following boundary consistency constraints must be satisfied:

[0041] S1.1 State Variable Initialization: In the initial stage of the multi-grid loop, the system initializes the local preprocessed displacement increment field of the current first-layer grid to zero in the memory of each GPU, and distributes the global residual passed from the outer conjugate gradient algorithm to the load residual field of the corresponding subdomain.

[0042] S1.2 High-frequency denoising of the physical field (smoothing iteration before activation): Activate the damped Jacobian smoother, set the expected number of iterations (typically 1 to 5) and the relaxation damping factor (e.g., 0.6). In each smoothing iteration, execute the following steps S1.3 to S1.5 sequentially: S1.3 Local Matrix-Free Template Operator Calculation (Solving for Current Internal Forces): To avoid explicit storage of a large stiffness matrix, the GPU initiates fine-grained parallel threads to perform matrix-free calculations. For each mesh element within the subdomain, the thread dynamically extracts the current displacement increments of its 24 degrees of freedom. Combining the Solid Isotropic Material Penalty Model (SIMP), the penalized stiffness scaling factor for the current element is calculated in real time using the element's current material density design variable, the preset minimum penalized elastic modulus, and the basic elastic modulus. Subsequently, this scaling factor is multiplied by the pre-calculated reference element stiffness matrix from the initialization phase, and a local multiply-add operation is performed with the mesh displacement vector (equivalent to obtaining the local internal force field).

[0043] S1.4 Asynchronous pipeline synchronization of boundary internal force fields (halo layer communication): When performing matrix-free computation in S1.3, a "boundary-first" asynchronous pipeline strategy is adopted to mask cross-GPU communication latency. (1) GPUs prioritize parallel computation of the mesh internal forces that are close to the physical partition boundary of the subdomain (i.e., the "halo layer").

[0044] (2) Once the boundary internal force calculation is completed, non-blocking PCIe bus communication is immediately triggered to send the internal force data of the halo layer to the adjacent left / right GPUs through the CPU memory buffer.

[0045] (3) While the data is being transmitted, the remaining computing power of the GPU continues to traverse and calculate the internal forces of the mesh within the subdomain (non-halo layer) in parallel.

[0046] (4) After the internal calculation is completed, the system performs communication barrier synchronization. The current GPU receives the internal force contribution of the aura layer from the neighboring GPU and safely accumulates it to the node of the physical boundary of the subdomain, thereby completing the seamless splicing of the global equivalent internal force field.

[0047] S1.5 Relaxation Update of Node Displacement Increments: After obtaining the complete global equivalent internal force field, the GPU restarts the parallel kernel to traverse all nodes in the current subdomain. The system reads the pre-assembled "inverse of the diagonal stiffness matrix" and uses the Jacobian update formula to correct the displacement increment of the current node: Latest displacement increment = Current displacement increment + Damping factor × Diagonal stiffness inverse matrix × (Load residual - Calculated global equivalent internal force). This completes one round of pre-smoothing, filtering out drastic geometric high-frequency errors in the mesh.

[0048] S1.6 Extracting the Downlink Residual Field: After reaching the preset number of pre-smoothing iterations, the physical field tends to be locally smooth. At this point, the GPU uses the final displacement increment field after multiple rounds of smoothing correction to perform one more complete S1.3 and S1.4 step (without matrix template multiplication and boundary aura synchronization) to obtain the smoothed final accurate internal force. Finally, each GPU subtracts this final accurate internal force from the load residual of its subdomain in parallel to obtain the denoised clean residual field. This residual field will serve as the external excitation input for the next layer (coarser mesh layer), preparing for the residual constraint and downlink operation in the subsequent step S2.

[0049] S2: The GPU uses a preset constraint operator to apply the first... The residuals of the first layer are constrained, and the resulting coarsened residuals are passed to the second layer. Layer meshes serve as load inputs for the next layer.

[0050] S2.1 3D Spatial Weighted Projection and Constraint Operations: The system initiates a dimensionality reduction mapping processing unit on the GPU side, constructing a dimensionality reduction mapping operation based on the geometric topological relationships of the 3D structured mesh, starting from the first... fine mesh up to the first Constraint operators for coarse mesh layers.

[0051] The GPU allocates computation threads to traverse the first... For each node in the coarse mesh layer, the system retrieves the neighboring nodes surrounding that node, centered on its spatial coordinates. A set of fine-mesh nodes. Based on the full-weighting principle or the inverse rule of trilinear interpolation, the system assigns spatial weights to the retrieved fine-mesh residual values ​​that are inversely proportional to the geometric distance. Through local multi-threaded parallel reduction, these weighted fine-mesh residual values ​​are accumulated onto the corresponding coarse-mesh nodes, thereby completing the spatial dimensionality reduction and feature aggregation of the residual vector within the physical subdomain.

[0052] S2.2 Coarse-grid scale aurora layer boundary synchronization: Due to the spatial span characteristics of the limiting operator, coarse-grid nodes adjacent to the GPU subdomain partitioning boundary cannot obtain the fine-grid residual contribution within the adjacent GPU subdomain during step S2.1, resulting in the boundary residual after local dimensionality reduction being in an "incompletely calculated" state. To ensure the global numerical conservation of the physics field, the system in the... At the coarse mesh scale of the layer, a cross-GPU boundary data collaboration mechanism is triggered: each GPU extracts the coarse mesh residual fragments that need to be passed outwards from the boundary of its own subdomain (halo layer). The system schedules the CPU's main memory buffer and performs non-blocking asynchronous boundary data transmission via a communication bus (such as PCIe), sending the residual fragments to the memory of adjacent GPU nodes. Each GPU receives the halo layer residual data from the adjacent subdomain and precisely superimposes it onto the residual value of the corresponding node on the physical boundary of its own subdomain. Through this cross-node overlapping summation, the system logically "stitches" together the cut physical boundaries, ensuring that the result of multi-GPU distributed constraint operations is absolutely equivalent to the result of single-GPU global computation.

[0053] S2.3 Hierarchical Depth Progression and State Handover: After completing the dimensionality reduction projection within the subdomain and the boundary data synchronization across subdomains, the system confirms the current subdomain's... The coarse-grained mesh has yielded a complete and accurate global equivalent residual field. Subsequently, the system updates the mesh level indicator variable of the control loop, making... This completes the state handover between the current coarse and fine mesh levels and prepares to transfer the computation flow to the next depth of multi-mesh recursive decision or preprocessing steps.

[0054] S3: Order .when Less than the preset bottom threshold When, perform pre-smoothing on the current layer and return to step S2 to continue the restriction operation; when equal When the coarsest mesh layer is reached, the direct solver is invoked to solve the [missing information - likely a specific problem or solution]. The state equation of the subdomain is obtained to obtain the first Exact solution for displacement of coarse mesh layer.

[0055] S3.1 Recursive condition determination and calculation flow: The system reads the currently updated grid level indicator variable. and compared with the system's preset coarsest grid level threshold. (i.e., the maximum depth of the bottom layer of the multi-grid V-loop) is compared and judged: like This indicates that the system is currently in the middle transition layer of a multi-grid system. At this point, the system is targeting the current... For each layer of mesh, repeat the pre-smoothing and matrix-free state equation calculation in step S1, and continue to execute step S2 to constrain and propagate the residuals to deeper levels.

[0056] like This indicates that the system has penetrated down to the coarsest grid layer (the bottom extremely coarse grid). At this scale, the number of grid nodes has been significantly reduced, and iterative smoothers are no longer effective in eliminating global low-frequency errors. The system then stops propagating downwards and enters the coarsest layer heterogeneous accurate solution stage, triggering the subsequent step S3.2.

[0057] S3.2 Data Upstream Convergence and Explicit Assembly of Sparse Matrices in the Coarsest Grid: In the coarsest layer (the... The system abandons the matrix-free computation mode on the GPU side and instead adopts an explicit matrix assembly strategy. Each GPU stops iterative smoothing operations, extracts the local subdomain residual vectors of the coarsest layer residing in its video memory, and sends them back to the CPU's host memory via a high-speed communication bus. The CPU seamlessly combines the local residual vectors scattered across multiple GPU nodes into a global right-hand side external load vector according to the topology numbering rules of the global mesh. Simultaneously, the system extracts the local stiffness properties of the coarsest layer mesh in the CPU memory and explicitly assembles them into a global ultra-coarse mesh stiffness matrix in Compressed Sparse Row (CSR) format.

[0058] S3.3 Direct and Accurate Solution Based on Preset Decomposition Factors: After completing the global load vector aggregation and stiffness matrix assembly, the system calls a high-precision sparse direct solver (such as the high-performance supernode Cholesky decomposition algorithm) on the CPU. To minimize the solution time of a single iteration, the system directly extracts the Cholesky decomposition factors (i.e., symbolic and numerical decomposition results) of the coarsest-layer stiffness matrix, which are pre-calculated and cached during the optimization initialization phase. Using these preset decomposition factors, the CPU performs forward and backward substitution solutions on the currently aggregated global load vectors. Through this direct solution process, the system obtains the globally accurate reference displacement correction solution at this level in one go, thus mathematically eliminating the residual persistent low-frequency rigid body displacement errors completely and accurately.

[0059] S3.4 Topology Decomposition and Memory Distribution of the Baseline Solution: After obtaining the globally accurate displacement correction solution on the CPU, the system, based on the initial physical subdomain segmentation mask, decomposes the globally continuous displacement vector into multiple locally corrected fragments corresponding one-to-one with the original GPU subdomains. Subsequently, the CPU distributes and copies these locally corrected fragments to the corresponding GPU's memory, completing the accurate update of the displacement state of the coarsest layer mesh. Simultaneously, the system initializes the error propagation ascending level indicator variable, setting $j = N$, to prepare for the subsequent "dimensionality-upgrading" calculation of upward interpolation.

[0060] S4: Order For the first The displacements calculated from the subdomains are mapped to the first subdomain using the prolongation operator. The layer undergoes error correction and post-smoothing processing; let The process continues until all subdomains have been traversed upwards, returning to layer 1 and obtaining the preprocessed residual vector. Then, the global residual of the layer 1 state equation is calculated (during which boundary data exchange between adjacent subdomains is triggered via OpenMP). If the global residual does not meet the preset convergence threshold, the search direction is updated and the process returns to step S2 for the next MGCG iteration; otherwise, the state equation solution converges, and the process proceeds to step S5.

[0061] S4.1 Up-interpolation and dimensionality-scaling of error correction: for the current... Layer (coarser mesh layer, initially) The system initiates the upward interpolation kernel in parallel on the GPU. Based on the three-dimensional trilinear interpolation operator (Prolongation Operator), the system will... The low-frequency displacement correction solution on the nodes of the coarse mesh layer is assigned spatial weights according to geometric topological distance and is equivalently superimposed onto the first layer. At the corresponding nodes of the layer (the finer first-level grid). Thus, the low-frequency mechanical response captured in the deep layer is accurately fed into the preprocessed displacement increment field of the shallower layer.

[0062] S4.2 Post-smoothing operation to filter spatial mapping noise: Due to the abrupt change in interlayer mesh resolution, geometric upscaling interpolation will inevitably cause noise in the first layer. Non-physical geometric high-frequency noise (numerical fluctuations) is excited in the displacement field of the layer. To correct this bias, after completing the interpolation mapping, the system immediately applies a correction to the first layer of each GPU. After smoothing iterations are initiated on the subdomains: each GPU again uses matrix-free template operators to calculate the local residuals of the current level; a damped Jacobi relaxation algorithm is used to filter high-frequency errors in the displacement increment field. During the post-smoothing process, the system synchronously triggers asynchronous pipelined communication across GPUs' Halo Layers to perform non-blocking exchange and overlap stitching of the smoothed physical boundary displacement data with adjacent subdomains, ensuring that the upgraded physical field quickly recovers smoothness and continuity at the cross-GPU boundaries.

[0063] S4.3 Upward Multi-Scale Backtracking and Preprocessing Direction Output: After completing the post-smoothing operation of the current layer, the system updates the upward level indicator variable. Determine the current level: If Then repeat steps S4.1 and S4.2 above, continuing to map and smooth the error correction amount layer by layer upwards; until... The system successfully backtracked to the finest original mesh layer. At this point, the single multigrid (MG) preprocessing loop was completely completed, and the system obtained the globally optimal conjugate gradient search direction after multi-scale full-band error correction.

[0064] S4.4 Global Convergence Verification of Outer-Layer State Equations: Based on the obtained optimal search direction and calculation step size, the system's main control thread globally updates the total displacement field of the 3D cantilever beam structure in the first layer (the finest mesh layer) and simultaneously updates the global external load residual vector. The system calculates the L2 norm of the global total residual in the current iteration step and compares it with the system's initial state residual to obtain the relative residual ratio: If this residual ratio is greater than the preset mechanical equilibrium tolerance threshold (e.g., ...), then... This indicates that the stress state of the current structure has not yet reached the finite element physical equilibrium. The system will return to step S1 with the latest residual field and start the next round of conjugate gradient and multi-grid preprocessing iterations. If the residual ratio is less than or equal to the preset mechanical equilibrium tolerance threshold, it is determined that the finite element state equation under the current topology has converged precisely, the internal physical field solution loop ends, and the optimization process continues to step S5.

[0065] S5: After the state equations converge, each GPU calculates in parallel the sensitivity of the material density of the first-layer subdomain mesh cells to the objective function (such as compliance). Based on the volume constraints and sensitivity information of the global mesh cells, the material density design variables of each mesh cell are updated using an optimization algorithm. After performing density filtering, the updated physical subdomain for the current optimization iteration step is obtained. The above steps are repeated until the entire topology optimization process reaches global convergence.

[0066] S5.1 Conjugate Gradient Iteration: Based on the V-cycle as a preconditioner, conjugate gradient iteration is performed. For each iteration, the preconditioning vector is first obtained through the V-cycle. Then calculate:

[0067] And update the search direction Intermediate vector Step length This leads to the updating of the displacement vector and residual vector until the condition is met.

[0068] S5.2 Compliance and Sensitivity Calculation: After the state equations converge, each GPU calculates compliance in parallel on its local subdomain.

[0069] And calculate the sensitivity of each unit:

[0070] The formula for calculating the sensitivity of the material density of the first-layer subdomain mesh element to compliance is as follows:

[0071] in, The overall compliance of the structure, As a penalty factor, The material density of the current unit. This represents the nodal displacement vector of the current element. The material density is The element foundation stiffness matrix at that time.

[0072] The CPU then retrieves the sensitivity and performs sensitivity filtering to suppress checkerboard and grid dependency issues.

[0073] S5.3 Design Variable Update: After sensitivity filtering is completed, each GPU performs design variable updates in parallel under volume constraints. An optimality criterion update strategy is preferred, using a binary search of Lagrange multipliers to obtain new design variables that meet the volume fraction requirements. To ensure iterative stability, a movement limit is set for each update. .

[0074] The formula for updating the material density of each mesh cell based on the volume constraint of the mesh cell and the sensitivity is as follows:

[0075] in, and These represent the relative material densities before and after the mesh cell update. The minimum density lower bound set for the singularity of the stiffness matrix. For the movement limit, The iteration factor is jointly determined by the objective function sensitivity and the volume constraint Lagrange multiplier. This is the numerical damping coefficient.

[0076] After the update, the CPU performs density filtering to obtain a new physical density field, which is used for solving the state equations in the next round.

[0077] S5.4 Convergence Criterion: Calculate the maximum change in design variables between two consecutive outer iterations:

[0078] when If the iteration stops and the topology is output, the iteration stops; otherwise, the iteration returns to the state equation solution step and continues. Large-scale topology optimization involves a large number of scalar reduction operations, such as right-hand side norm, residual norm, rho, and total volume calculation. This paper implements a "device-side local reduction + host-side aggregation" approach to handle these scalars: each GPU first performs parallel reduction on its local subdomain, and then aggregates the local results into a global scalar through atomic updates or OpenMP reduction. For scalars that are more sensitive in long iteration chains, the code provides a mixed-precision implementation, using double precision for host-side accumulation to enhance iteration stability. During the design variable update phase, the program performs volume fraction checks, moving limit constraints, and design variable updates under binary search in parallel on the GPU side, and then sends the intermediate design variables back to the CPU side for density filtering. Thus, although the OpenMP framework has migrated the main high-density numerical operations to the device side, it still retains the filtering and control links involving the CPU. This choice helps control implementation complexity, but it also introduces an additional data backhaul path.

[0079] like Figure 2 As shown, this illustrates the spatial (across GPUs) and scale (across layers) two-dimensional data flow relationship of the Multigrid algorithm under a dual-GPU architecture.

[0080] Vertical direction (scale transformation): The blue "Restriction" arrow indicates that the residual is transferred from the fine mesh to the coarse mesh; the green "Prolongation" arrow indicates that the correction amount is mapped from the coarse mesh to the fine mesh.

[0081] Horizontal direction (spatial synchronization): A double-headed arrow running through the dashed line (physical partition boundary) indicates that "smoothing and halo exchange" is being performed between adjacent GPUs at both the fine and coarse mesh levels.

[0082] This diagram illustrates that the algorithm not only spatially segments the model (assigning it to left / right subdomains) but also rigorously performs communication of aura layer data at every resolution level of the multi-mesh. This ensures that regardless of mesh coarsening or refinement, cross-GPU computation results are always absolutely equivalent to single-GPU computation results, guaranteeing the continuity and conservation of the global physics field.

[0083] like Figure 3As shown, the complete control flow and hardware / software mapping of the solver in a single V loop iteration are depicted. The computational flows on both sides (GPU 0 and GPU 1) are symmetrically parallel. The diagram illustrates the downward process from "pre-smoothing" → "halo swap" → "constraint operator (dimensionality reduction)" and the upward process of "interpolation operator" $\rightarrow$ "post-smoothing". At the bottom of the diagram, when it is determined that the coarsest layer has been reached ( When the GPU computing flows on both sides converge to the center, the CPU performs "coarse mesh convergence and direct solution".

[0084] A vertically aligned domain decomposition architecture was established. For the global computation domain, the physical space partitioning boundary remains consistent regardless of the resolution. Assuming the global domain is divided into two subdomains along the X-axis, if a physical point is within the subdomain managed by the first GPU, all corresponding node data from the finest to the coarsest grid at that location are allocated to the GPU's video memory. As shown in the figure, this design ensures that the physical spatial hierarchy does not change with the grid resolution. Thanks to this, the computation process of constraint and interpolation operators is completely transformed into on-device video memory operations: when performing constraint operations, the GPU only needs to read the subset of fine grids in its local video memory to complete the aggregation of coarse grid residuals; similarly, the interpolation operation only needs to read the local coarse grid correction amount to update the fine grid solution. The entire V-Cycle uplink and downlink process is completed in parallel within the GPU, without triggering expensive PCIe bus cross-device data migration, and cross-GPU MPI communication is strictly limited to smooth operations at the same level.

[0085] like Figure 4 As shown, the same 3D cantilever beam design problem is illustrated at four different mesh resolutions ( , , , The final optimized topology array is compared. Each row represents a resolution, and the columns show different viewing angles (front view, side view, top view, axonometric).

[0086] The figure illustrates the framework's ability to solve 3D topology optimization problems of varying scales.

[0087] like Figure 5 As shown, High-resolution magnified close-up of the local optimization results (scale, i.e., more than 33 million grid cells and hundreds of millions of degrees of freedom).

[0088] This demonstrates that, under such a vast degree of freedom, the optimized 3D cantilever beam exhibits an extremely smooth surface, free from any numerical instability common in topology optimization. This proves that the multi-GPU heterogeneous parallel architecture proposed in this patent truly possesses the capability to solve industrial-grade, multi-hundred-million-degree-of-freedom ultra-large-scale topology optimization problems on a single-node machine.

[0089] Those skilled in the art will readily understand that the above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A CPU / multi-GPU heterogeneous parallel acceleration method for topology optimization of three-dimensional cantilever beams based on OpenMP, characterized in that, The method includes the following steps: The subdivided mesh design domain of the 3D cantilever beam to be optimized is divided into multiple subdomains. Topology optimization is performed on each subdomain. The optimized subdomains are then joined together to obtain the optimized 3D cantilever beam design domain. The optimization of each subdomain is carried out according to the following steps: S1 i=1, take the grid cells of the current subdomain as the i-th layer, and calculate the residuals of each grid cell in the current subdomain; S2 updates the residuals of each grid cell in the current subdomain by multiplying the preset constraint operator with the residual, and updates each grid cell according to the updated residuals, thereby updating the current subdomain; S3 i=1+1, when i is less than or equal to the preset threshold N, return to step S2; otherwise, solve the residual of each grid cell in the Nth subdomain; S4 j=N, use the residuals of each grid cell in the j-th subdomain to reprocess the subdomain. This reprocessing includes interpolation, error correction and post-smoothing. j=j-1, until the reprocessing of all subdomains is completed. Calculate the residuals of the state equations of the grid cells in the first subdomain. If the residuals do not meet the preset threshold, return to step S2. Otherwise, proceed to S5. S5 Calculate the sensitivity of the material density of the first-layer subdomain mesh cell to compliance; update the material density of each mesh cell according to the volume constraint of the mesh cell and the sensitivity, and use the updated material density to update the first-layer subdomain mesh cell. The updated first-layer subdomain mesh cell is the optimized subdomain.

2. The method for topology optimization of three-dimensional cantilever beams based on OpenMP and heterogeneous parallel acceleration using CPU / multi-GPU architectures as described in claim 1, characterized in that... After dividing the subdivided mesh design domain of the three-dimensional cantilever beam to be optimized into multiple subdomains, each subdomain is divided into an aura layer and a non-aura layer. Then, each subdomain is mapped to the GPU. Next, topology optimization is performed on each subdomain in the GPU. In step S3, when solving the state equation of the Nth layer subdomain, the design variables of the Nth layer need to be passed to the CPU connected to the GPU. The residuals of each mesh cell of the Nth layer subdomain are solved in the CPU.

3. The method for topology optimization of three-dimensional cantilever beams based on OpenMP and heterogeneous parallel acceleration using CPU / multi-GPU architectures, as described in claim 2, is characterized in that... In step S1, the residuals of each grid cell are obtained by first solving the state equation of the subdomain to obtain the displacement of each grid cell, and then using the displacement to solve for the residuals.

4. The method for topology optimization of three-dimensional cantilever beams based on OpenMP and heterogeneous parallel acceleration using CPU / multi-GPU architecture, as described in claim 3, is characterized in that... The state equation of the subdomain is as follows: KU=F in, The subdomain stiffness matrix is ​​related to the material distribution. and These are the subdomain displacement vector and the subdomain external load vector, respectively.

5. The method for topology optimization of three-dimensional cantilever beams based on OpenMP and heterogeneous parallel acceleration using CPU / multi-GPU architectures, as described in claim 4, is characterized in that... The formula for calculating the residual is as follows: in, It represents the residual.

6. A method for topology optimization of three-dimensional cantilever beams based on OpenMP and heterogeneous parallel acceleration using CPU / multi-GPU architectures, as described in claim 3 or 4, characterized in that... When calculating the state equation of a subdomain, for subdomains from layer 1 to N-1, the data of the aura layer in adjacent subdomains need to be exchanged in the CPU, and then the state equation of the current subdomain is calculated using the data of the aura layer in the adjacent subdomains and the non-aura layer in the current subdomain. For the Nth subdomain, the data of the aura layer and the non-aura layer in the current subdomain need to be transferred to the CPU for joint calculation to obtain the state equation of the current subdomain.

7. The method for topology optimization of three-dimensional cantilever beams based on OpenMP and heterogeneous parallel acceleration using CPU / multi-GPU architectures as described in claim 1, characterized in that... In step S4, the interpolation, error correction, and post-smoothing processes employ trilinear interpolation, direct summation, and damped Jacobian algorithms, respectively.

8. The method for topology optimization of three-dimensional cantilever beams based on OpenMP and heterogeneous parallel acceleration of CPU / multi-GPU as described in claim 1, characterized in that, In step S5, the formula for updating the material density of each mesh cell based on the volume constraint of the mesh cell and the sensitivity is as follows: in, and These represent the relative material densities before and after the mesh cell update. The minimum density lower bound set for the singularity of the stiffness matrix. For the movement limit, The iteration factor is jointly determined by the objective function sensitivity and the volume constraint Lagrange multiplier. This is the numerical damping coefficient.

9. The method for topology optimization of three-dimensional cantilever beams based on OpenMP and heterogeneous parallel acceleration using CPU / multi-GPU architecture as described in claim 1, characterized in that, In step S5, the formula for calculating the sensitivity of the material density of the first-layer subdomain mesh element to compliance is as follows: in, The overall compliance of the structure, As a penalty factor, The material density of the current unit. This represents the nodal displacement vector of the current element. The material density is The element foundation stiffness matrix at that time.

10. A CPU / multi-GPU heterogeneous parallel acceleration system for three-dimensional cantilever beam topology optimization based on OpenMP, characterized in that, The system includes an actuator for performing a CPU / multi-GPU heterogeneous parallel acceleration method for three-dimensional cantilever beam topology optimization based on any one of claims 1-9.

Images