PD-FEM coupling implicit solution method based on GPU parallel

Abstract

The invention is suitable for the technical field of computational mechanics and high-performance computation crossing, and relates to a PD-FEM coupling implicit solution method based on GPU parallelism, which comprises the following steps: S10, model preprocessing and data initialization; s20, assembling an integral stiffness matrix based on cooperation of the CPU and the GPU; s30, carrying out parallel solution on the system linear equation set based on the GPU; s40, carrying out GPU-based damage evolution parallel calculation; s50, damage iteration control and stiffness matrix parallel updating are carried out; and S60, result return and post-processing are carried out. The method is simple in process and convenient to operate, not only obtains breakthrough progress in calculation speed, but also shows remarkable advantages in simulation precision, numerical stability and application range, and provides an efficient, reliable and powerful numerical simulation scheme for solving the problem of complex structure fracture analysis existing in the engineering field for a long time.

Description

Technical Field

[0001] This invention belongs to the field of computational mechanics and high-performance computing, and in particular relates to a GPU-parallel PD-FEM coupled implicit solution method. Background Technology

[0002] Accurate simulation of crack propagation is crucial for the safety assessment and life prediction of engineering structures. Currently, mainstream numerical methods include the Finite Element Method (FEM) and Perifield Dynamics (PD). The FEM, a mature method based on continuum mechanics, demonstrates high efficiency and accuracy in handling continuous problems such as elastic deformation and has been widely applied in engineering. However, because its theoretical foundation relies on the spatial partial derivatives of the displacement field, singularity problems arise when dealing with discontinuities such as crack tips. External fracture criteria and complex mesh reconstruction techniques are required to simulate crack propagation, which not only increases algorithm complexity but may also introduce human error, affecting the accuracy of crack path prediction. In contrast, PD is an integral theory based on nonlocality. It describes interaction forces through "bonds" between point pairs, avoiding dependence on spatial derivatives. Therefore, it can naturally and spontaneously simulate crack initiation and propagation without the need for complex fracture criteria. However, the nonlocality of the near-field dynamics method requires high-density discretization of the model, which results in a computational cost far exceeding that of the traditional finite element method. When applied to the overall analysis of large structures, its computational cost is very high, and its computational efficiency is affected.

[0003] Patent application CN119864107A discloses a method for predicting crack propagation in elastobrittle materials based on bonded peri-field dynamics, comprising the following steps: S1, establishing a geometric model of the elastobrittle material structure and setting boundary conditions as a preprocessing input file; S2, establishing equilibrium equations on the elements based on the finite element method; S3, solving the equilibrium equations and calculating the damage to the elastobrittle material matrix; S4, displaying the calculation results through a post-processing platform. This patent application uses bonded peri-field dynamics to simulate the entire process of crack initiation, propagation, and convergence in elastobrittle materials. However, it involves a large computational load and is inefficient, exhibiting the same drawbacks as existing technologies.

[0004] Therefore, how to provide a numerical analysis method that combines computational accuracy, numerical stability, and computational efficiency is a problem that urgently needs to be solved by researchers in this field. Summary of the Invention

[0005] To address the shortcomings of existing technologies, the purpose of this invention is to provide a GPU-parallel PD-FEM coupled implicit solution method to solve the problems of low computational accuracy and efficiency in existing numerical analysis methods.

[0006] To solve the above-mentioned technical problems, the present invention adopts the following technical solution:

[0007] This invention provides a GPU-parallel PD-FEM coupled implicit solution method, comprising the following steps:

[0008] S10. Define the geometric and material properties of the model, discretize it and divide it into PD and FEM sub-regions, construct the bond connection relationship between PD material points, and transfer all the basic data structures required for calculation from the host memory to the GPU device memory.

[0009] S20. Construct the overall system stiffness matrix K in the initial state. The stiffness contributions of the FEM unit and the interface coupling unit with smaller computational load are calculated on the CPU, while the stiffness contributions of all bonds in the PD region with huge computational load are calculated and assembled in large-scale parallel by starting the GPU kernel function, and finally forming a complete stiffness matrix in the GPU device memory.

[0010] S30. Use a parallelized iterative solver to solve the large sparse linear equation system Ku=F stored in the GPU device memory to obtain the displacement vector u of all FE nodes and PD material points in the current system state.

[0011] S40. By starting the GPU kernel function, the elongation of each bond is calculated in parallel and compared with the material's fracture criterion. If the fracture condition is met, the bond is marked as fractured in the global bond state array in the GPU device's memory.

[0012] S50: The host CPU checks whether any new damage has occurred in S40. If so, the system stiffness changes. At this time, the GPU kernel function is started, and the stiffness matrix K is reassembled in parallel according to the information in the key state array. After the update is completed, the process returns to step S30 and solves the problem again using the new stiffness matrix. If no new damage has occurred, the convergence of the residual force of the current iteration is further judged. If the residual force meets the preset threshold, the system is considered to have reached equilibrium and the iteration ends. Otherwise, the process returns to step S30 to continue iterating and solving the problem.

[0013] S60. After the damage iteration cycle ends, the final calculation results, including the displacement field and damage distribution, are transferred from the GPU device memory back to the host memory for subsequent visualization analysis and data archiving.

[0014] Furthermore, the specific steps of S10 are as follows:

[0015] S101. Define the geometric model, material properties, boundary constraints, and external load conditions of the analysis object. Discretize the geometric model spatially, dividing the model into FEM sub-regions and PD sub-regions. In the FEM sub-region, generate a standard finite element mesh containing nodes and elements. In the PD sub-region, generate discrete material points.

[0016] S102. Define an interface region at the junction of the FEM subregion and the PD subregion. Within this region, establish the coupling connection between the FEM node and the PD material point.

[0017] S103. Traverse each material point i in the PD sub-region and search for all other material points in its near-field neighborhood. Store all the searched material point pairs as a global key list.

[0018] S104. In the host memory, based on the results of discretization and neighborhood search, construct all necessary data arrays, including the coordinate arrays of nodes and material points, the FEM element connection relationship array, the interface coupling element connection relationship array, and the PD bond list generated in S103. In the GPU device memory, allocate contiguous memory space for parallel computing to store the overall stiffness matrix K, displacement vector u, load vector F, global bond state array, material parameters, and coordinates.

[0019] S105. Perform a one-time data copy operation from host memory to GPU device memory, transferring the coordinate array, all connection relationship arrays, material parameters, and initialized load vector F constructed in S104 from the CPU to the corresponding memory space of the GPU.

[0020] Furthermore, S102 specifically includes:

[0021] S121. In the discretized model, explicitly define the FEM sub-region and the PD sub-region, identify the interface located at the junction of the two sub-regions, and determine the FE nodes and PD material points on the interface.

[0022] S122. Introduce truss elements to bridge the approach field dynamics sub-region and finite element sub-region;

[0023] S123. Calculate the element stiffness matrix of each truss element using the standard finite element method. .

[0024] Furthermore, the specific steps of S20 are as follows:

[0025] S201. On the CPU, calculate the stiffness contribution of the FEM sub-region and the interface region:

[0026] FEM subregion stiffness: The CPU iterates through all FEM elements and calculates the element stiffness matrix K for each element. fem ;

[0027] Interface coupling stiffness: The CPU traverses all interface truss elements established in S122 and calculates the global stiffness matrix K of each coupling element according to the method in S123. t ;

[0028] Calculate all K fem and K t The index position information of the matrix and its position in the overall stiffness matrix K is stored in a temporary buffer in the CPU memory, ready for subsequent transmission to the GPU;

[0029] S202. On the GPU device, the stiffness contribution of the PD sub-region is calculated and assembled in parallel, which is achieved by starting the GPU kernel function, and a GPU thread is allocated for each key in the key list in S103.

[0030] S203, The stiffness contribution data K calculated by the CPU in S201 is... fem and K t The stiffness contribution from the CPU is transferred to the GPU device memory, and a GPU kernel function is started. The stiffness contribution from the CPU is then added to the overall stiffness matrix K in the GPU device memory through atomic operations, thus completing the final assembly.

[0031] Furthermore, the specific steps of S30 are as follows:

[0032] S301. Initialize the vectors required for solving the problem on the GPU, including the initial displacement vector u0 and the residual vector r0 = F - K. u0 And the initial search direction vector p0 = r0;

[0033] S302. Execute an iterative loop on the GPU until the convergence condition is met;

[0034] S303. After the iteration converges, the final global displacement vector u is stored in the GPU device memory.

[0035] Furthermore, in S302, each iteration includes the following GPU parallel computation:

[0036] Sparse matrix-vector multiplication: Parallel computation This operation utilizes sparse matrices and vectors p stored in GPU memory. k Multiplication and accumulation operations are performed concurrently by a large number of GPU cores;

[0037] Parallel reduction: The inner product is calculated using a GPU-based parallel reduction algorithm to determine the step size. ;

[0038] Vector update: The displacement solution vector and residual vector are updated in parallel, with each thread responsible for one or more components. and ;

[0039] Convergence criterion: The norm of the new residual is calculated in parallel, and then calculated through another parallel reduction. The iteration is then compared with a preset convergence tolerance; if the result is less than the tolerance, the iteration stops.

[0040] Direction vector update: If convergence is not achieved, the search direction for the next iteration is calculated and updated through parallel reduction and AXPY operations. ,in .

[0041] Furthermore, the specific steps of S40 are as follows:

[0042] S401, The host CPU determines the total length N of the key list in S103. bonds Start a configuration with N bonds Each GPU thread is uniquely mapped to a PD key;

[0043] S402, On the GPU, all N bonds Several threads execute concurrently. Each thread k first reads the displacement u of the corresponding material points i and j at both ends of the bond from the global displacement vector u in the GPU device memory. i and u j Subsequently, the thread determines the initial relative position of the material points. and current relative displacement Calculate the relative elongation of the bond. ;

[0044] S403. Thread k compares its calculated elongation s with the material's critical elongation s0. If s > s0 and the bond is still in an intact state in the global bond state array, then the thread updates the bond's state to broken.

[0045] Furthermore, the specific steps of S50 are as follows:

[0046] S501. The GPU performs a parallel reduction operation to quickly count the number of keys whose states have changed in the global key state array, and obtains a new damage count value. The host CPU retrieves the new damage count value from the GPU.

[0047] S502, The host CPU judges the newly added damage count:

[0048] If the new damage count > 0, it indicates that the stiffness of the system has changed due to material damage, and the current displacement solution u is no longer valid. The system must update the stiffness matrix and solve it again, and execute step S503.

[0049] If the newly added damage count is 0, and after further convergence judgment on the residual force of the current iteration, the residual force still meets the preset threshold, it indicates that under the current load, the material damage no longer expands, the system has reached equilibrium, the implicit iteration loop ends, and step S60 is executed.

[0050] S503, Parallel update of stiffness matrix to reconstruct the overall system stiffness matrix K;

[0051] S504. After the GPU completes the update of the stiffness matrix K in S503, the program control flow returns to S30. Based on the GPU-based parallel solution step for the system's linear equations, the system uses the updated K matrix, combined with the original load vector F, to re-solve Ku=F, thereby obtaining the new displacement field u. new And return to S40.

[0052] Furthermore, S503 specifically includes:

[0053] S531. Before reassembly, clear all stiffness contributions related to PD bonds in the overall stiffness matrix K stored in the GPU device memory.

[0054] S532. The host CPU starts a GPU kernel function and allocates a GPU thread for each key in the key list. In this kernel function, each thread k corresponds to key (i, j) and executes a conditional assembly logic:

[0055] Reading the state: Thread k first reads its state value from the global key state array in the GPU global memory;

[0056] Conditional judgment: If the state is intact, thread k calculates the global stiffness matrix k of the key. bond And use the atomicAdd() function to add its 36 components to the corresponding positions of the stiffness matrix K;

[0057] If the state is fractured, thread k will skip all subsequent stiffness calculations and atomicAdd() assembly operations and terminate the task directly. After this process is completed, the K matrix in the GPU device memory will be updated to a new stiffness matrix that contains only all intact bonds and contributions from FEM units and interface units.

[0058] Furthermore, the specific steps of S60 are as follows:

[0059] S601, The host CPU initiates a data transfer operation to copy the key result data stored in the GPU device memory back to the CPU host memory. The data includes the final global displacement vector u and the final global key state array.

[0060] S602: After receiving the data, the host CPU performs post-processing operations, including saving the displacement and damage data to disk files, and calling the visualization program library to graphically display the displacement field and global key state array. It can also extract key data points for quantitative analysis as needed.

[0061] S603. After post-processing is completed, the entire simulation process ends.

[0062] The GPU-parallel PD-FEM coupled implicit solution method provided by this invention has at least the following advantages compared with existing technologies:

[0063] Existing numerical analysis methods, including the finite element method (FEM) and peridynamics (PD), suffer from low computational accuracy and efficiency. This invention offers a simple and convenient process, aiming to accurately and efficiently simulate the entire process of solid materials under external loads, from elastic deformation to crack initiation, propagation, and eventual failure. It applies peridynamics (PD) theory to the predicted potential fracture region and the finite element method (FEM) to the remaining linear elastic region, and designs a stable and efficient interface coupling mechanism. Its core innovation lies in addressing the high computational overhead caused by damage evolution in the implicit solution format of the coupled model by proposing a complete GPU parallel computing scheme. This migrates computationally intensive tasks from the CPU to the GPU, fundamentally solving the problem of balancing accuracy and efficiency in fracture simulation in existing technologies, and significantly improving the computational accuracy and efficiency of numerical analysis methods compared to existing techniques. This invention not only achieves a breakthrough in computing speed, but also demonstrates significant advantages in simulation accuracy, numerical stability, and application scope. It provides an efficient, reliable, and powerful numerical simulation solution for solving the long-standing problem of fracture analysis of complex structures in the engineering field. Attached Figure Description

[0064] To more clearly illustrate the solution of the present invention, a brief introduction will be given to the drawings used in the description of the embodiments below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0065] Figure 1A flowchart of a GPU-parallel PD-FEM coupled implicit solution method provided in an embodiment of the present invention;

[0066] Figure 2 A schematic diagram of the PD-FEM interface coupling mechanism of a GPU-parallel PD-FEM coupled implicit solution method provided in an embodiment of the present invention;

[0067] Figure 3 A schematic diagram illustrating the parallel assembly principle of the stiffness matrix in a GPU-parallel PD-FEM coupled implicit solution method provided in an embodiment of the present invention;

[0068] Figure 4 A CPU-GPU collaborative computing task allocation diagram for the implicit iterative process in a GPU-parallel PD-FEM coupled implicit solution method provided in an embodiment of the present invention;

[0069] Figure 5 A schematic diagram of the internal structure of an electronic device capable of executing a GPU-parallel PD-FEM coupled implicit solution method, provided for embodiments of the present invention;

[0070] Figure 6 The initial model setup diagram for a GPU-parallel PD-FEM coupled implicit solution method provided in an embodiment of the present invention;

[0071] Figure 7 A schematic diagram of crack propagation in a GPU-parallel PD-FEM coupled implicit solution method provided in an embodiment of the present invention;

[0072] Figure 8 This diagram illustrates the comparison between the load-CMOD curve and the theoretical solution in a GPU-parallel PD-FEM coupled implicit solution method provided in this embodiment of the invention. Detailed Implementation

[0073] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art. The terminology used in this specification is for the purpose of describing particular embodiments only and is not intended to limit the invention. For example, terms such as “length,” “width,” “upper,” “lower,” “left,” “right,” “front,” “rear,” “vertical,” “horizontal,” “top,” “bottom,” “inner,” and “outer” indicate orientations or positions based on the accompanying drawings and are for ease of description only, and should not be construed as limiting the technical solution. The terms “comprising” and “having,” and any variations thereof, in the specification, claims, and accompanying drawings of this invention, are intended to cover non-exclusive inclusion. The terms “first,” “second,” etc., in the specification, claims, and accompanying drawings of this invention are used to distinguish different objects and not to describe a particular order. In the specification, claims, and accompanying drawings of this invention, when an element is referred to as “fixed to,” “mounted to,” “set on,” or “connected to” another element, it may be directly or indirectly located on that other element. For example, when an element is referred to as “connected to” another element, it may be directly or indirectly connected to that other element. Furthermore, the reference to "embodiment" herein means that a particular feature, structure, or characteristic described in connection with an embodiment may be included in at least one embodiment of the invention. The appearance of this phrase in various places throughout the specification does not necessarily refer to the same embodiment, nor is it a separate or alternative embodiment mutually exclusive with other embodiments. It will be explicitly and implicitly understood by those skilled in the art that the embodiments described herein can be combined with other embodiments.

[0074] This invention provides a GPU-parallel PD-FEM coupled implicit solution method, applied to the safety assessment and life prediction of engineering structures. The GPU-parallel PD-FEM coupled implicit solution method is as follows:

[0075] S10. Define the geometric and material properties of the model, discretize it, and divide it into PD and FEM sub-regions. Construct the bond connections between PD material points and transfer all the basic data structures required for computation from host memory to GPU memory. S20. Construct the overall system stiffness matrix K in the initial state. The stiffness contributions of the computationally less computationally intensive FEM units and interface coupling units are calculated on the CPU, while the stiffness contributions of all bonds in the computationally intensive PD region are calculated and assembled in large-scale parallel computation using GPU kernel functions, ultimately forming a complete stiffness matrix in GPU memory. S30. Use a parallelized iterative solver to solve the large sparse linear equation system Ku=F stored in GPU memory to obtain all F values ​​in the current system state. The displacement vector u of node E and material point PD; S40, by starting the GPU kernel function, the elongation of each bond is calculated in parallel and compared with the fracture criterion of the material. If the fracture condition is met, the bond is marked as fractured in the global bond state array in the GPU device memory; S50, the host CPU checks whether new damage has occurred in S40. If so, the system stiffness changes. At this time, the GPU kernel function is started, and the stiffness matrix K is reassembled in parallel according to the information in the bond state array. After the update is completed, it returns to step S30 and solves again using the new stiffness matrix. If no new damage has occurred, the convergence of the residual force of the current iteration is further judged. If the residual force meets the preset threshold, the system is considered to have reached equilibrium and the iteration ends. Otherwise, it returns to step S30 to continue iterating; S60, after the damage iteration cycle ends, the final calculation results, including the displacement field and damage distribution, are transferred from the GPU device memory back to the host memory for subsequent visualization analysis and data archiving.

[0076] This invention features a simple process and convenient operation. It has achieved a breakthrough in calculation speed, and also shows significant advantages in simulation accuracy, numerical stability and application scope. It provides an efficient, reliable and powerful numerical simulation solution for solving the long-standing problem of fracture analysis of complex structures in the engineering field.

[0077] To enable those skilled in the art to better understand the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings.

[0078] This invention provides a GPU-parallel PD-FEM coupled implicit solution method, applied to the safety assessment and life prediction of engineering structures, combined with... Figures 1 to 8 In this embodiment, the GPU-parallel PD-FEM coupled implicit solution method is as follows:

[0079] S10. Model Preprocessing and Data Initialization: This stage is completed on the host (CPU), including defining the geometric and material properties of the model, discretizing it and dividing it into PD and FEM sub-regions, constructing the "bond" connection relationship between PD material points, and finally transferring all the basic data structures required for calculation from the host memory to the GPU device memory.

[0080] Specifically, in this embodiment, step S10 is a model preprocessing and data initialization method, designed to prepare all necessary geometric information, material parameters, topological connections, and GPU memory data for subsequent CPU-GPU collaborative computing. This method is executed on the host (CPU), and the specific steps are as follows:

[0081] S101. Model Definition and Discretization: First, define the geometric model of the analysis object, material properties (including FEM parameters such as elastic modulus and Poisson's ratio, and PD parameters such as near-field dynamic critical elongation), boundary constraints, and external load conditions. Second, spatially discretize the geometric model. Based on the analysis requirements, divide the model into a finite element (FEM) sub-region where continuity and linear elasticity are expected, and a near-field dynamic (PD) sub-region where damage and fracture are expected. Within the FEM sub-region, generate a standard finite element mesh containing nodes and elements (such as quadrilateral or hexahedral elements). Within the PD sub-region, generate discrete material points.

[0082] S102. Interface Region Definition and Coupling: At the boundary between the FEM sub-region and the PD sub-region, an interface region is defined. Within this region, the coupling connection between the FEM node and the PD material point is established. Specifically, the coupling connection between the interface FEM node and its near-field neighborhood is defined. Truss elements are established between PD material points within the range. These truss elements will act as bridges for force and displacement transmission, ensuring deformation coordination between the two sub-regions.

[0083] S103, PD Near-Field Neighborhood Search and Key List Construction: This step is the core of constructing the topological relationships of the PD model. It aims to identify all interaction pairs (i.e., "bonds") between PD matter points, traversing each matter point i in the PD sub-region and searching its near-field neighborhood (a neighborhood with a radius of...). All other material points within the spherical region. All searched material point pairs are stored as a global "Bond List," which serves as the foundational data structure for subsequent GPU parallel computing. This data structure is optimized for large-scale parallel access by GPU threads.

[0084] S104, GPU Data Structure Initialization and Memory Allocation: In the host (CPU) memory, based on the discretization and neighborhood search results, all necessary data arrays are constructed, including: coordinate arrays of nodes and material points, arrays of FEM element connections, arrays of interface coupling element connections, and the PD "bond list" generated in the previous step. In the GPU device memory, large blocks of contiguous memory space are allocated for the parallel computations to be performed. These spaces are used to store: the global stiffness matrix K: typically allocated using a sparse matrix storage format suitable for parallel computation; displacement vector u and load vector F: allocating space for all degrees of freedom in the system (including FE nodes and PD material points); the global bond status array: this is a crucial array whose length is equal to the total number of bonds in the PD "bond list". This array is initialized to an "intact" state (all element values ​​are 1) to mark whether each bond is broken in step S40; and other auxiliary arrays required for computation, such as material parameters and coordinates.

[0085] In this embodiment, the overall stiffness matrix K stored in the GPU device memory can adopt the Compressed Sparse Row (CSR) format. In other embodiments, to optimize the execution efficiency of the sparse matrix-vector multiplication (SpMV) in step S30 on the GPU, other GPU-friendly sparse storage formats can be selected based on the sparsity characteristics of the matrix. These formats include: Coordinate (COO) format, ELLPACK (ELL) format, or a hybrid format that combines the advantages of ELL and COO.

[0086] S105 Initial Data Transfer: Perform a one-time data copy operation from the host memory to the GPU device memory. The coordinate array, all connection relationship arrays (FEM elements, interface elements, PD key list), material parameters, and the initialized load vector F constructed in S104 are transferred from the CPU to the corresponding memory space on the GPU. After this step is completed, the GPU has all the initial data required to execute S20 and subsequent calculation steps, thereby minimizing data communication between the CPU and GPU in subsequent iterative calculations.

[0087] Furthermore, in this embodiment, step S102 is an interface region and coupling method, which aims to establish a stable transition region that can accurately transmit force and displacement, so as to ensure deformation coordination and force balance between the continuous medium field of the FEM sub-region and the discrete material point field of the PD sub-region. Its core idea is to "bridge" the two different types of degrees of freedom by introducing virtual interface units. The method specifically includes the following steps:

[0088] S121. Interface layer identification: In the discretized model, the FEM sub-region and PD sub-region are clearly defined, the "interface" located at the junction of the two sub-regions is identified, and the "FE node" and "PD material point" on the interface are determined.

[0089] S122. Establish coupling connection: Introduce truss elements to bridge the approach field dynamics sub-region and the finite element sub-region, such as... Figure 2 As shown, at the interface, an FE node will be related to its scope (as shown in the figure). Multiple PD material points within the circle shown are connected. This connection is achieved by creating "bar elements" (i.e., truss elements), with each truss element connecting an FE node and a PD material point.

[0090] S123. Coupled Element Stiffness Calibration: The introduced truss element is treated as a special type of finite element. The element stiffness matrix of each truss element is calculated according to the standard finite element method. When assembling the overall stiffness matrix K in step S20, the stiffness matrices of these truss elements are... Through coordinate transformation and assembly, the coordinates are superimposed onto the corresponding degrees of freedom positions of the system's total stiffness matrix K. Through this coupling method, an explicit stiffness connection is established between the degrees of freedom of the FEM nodes and the degrees of freedom of the PD material points, realizing the mechanical coupling of the two sub-regions.

[0091] Furthermore, in this embodiment, in the coupling element stiffness calibration of step S123, to ensure the correct transfer of force and displacement between the FEM sub-region and the PD sub-region, this calculation method aims to correctly convert the axial stiffness of the truss element and contribute it to the overall system stiffness matrix in the global coordinate system. The specific implementation steps are as follows:

[0092] a. Establishing a local coordinate system: Consider a single truss element connecting node i (FE) and material point j (PD). Establish a local coordinate system x', y', z'. Make the x' axis of this local coordinate system align with the axial direction of the element, pointing from point i to material point j.

[0093] b. Definition of Local Stiffness Matrix: In this local coordinate system x'y'z', the mechanical behavior of the truss element is simplified, and it only has stiffness in the axial direction. The displacement vectors of nodes i and j in the local coordinate system can be expressed as u' = [u i ', v i ',w i ', u j ', v j ', w j '] T The corresponding local force vector is f' = [f x'i , fy'i , f z'i , f x'j , f y'j , f z'j ] T The local stiffness matrix K of the truss element p It has non-zero terms only in the degrees of freedom related to x, and can be expressed as:

[0094] Among them, K p It is the axial stiffness scalar of the truss unit.

[0095] c. Direction Cosine Calculation: To transform the local stiffness matrix to the global coordinate system x, y, z, it is necessary to calculate the direction cosines l, m, n of the local x' axis in the global coordinate system. Let the global coordinates of node i be (x... i y i , z i The global coordinates of matter point j are (x... j y j , z j ), the element length is .

[0096] The formula for calculating the direction cosine is as follows:

[0097] .

[0098] d. Global stiffness matrix construction: Utilizing standard finite element coordinate transformation theory, based on direction cosines l, m, n and axial stiffness K... p The truss element can be directly constructed in the global coordinate system x, y, z. The element stiffness matrix k is 6. This matrix k establishes the global displacement vector [u i , v i , w i , u j , v j , w j ] T The relationship with the global force vector.

[0099] The specific form of matrix k is:

[0100]

[0101] e. Assemble to the overall system stiffness matrix: The global stiffness matrix k of each interface truss element calculated in step d is added to the overall system stiffness matrix K described in step S2 according to the index positions of the FE node i and PD material point j connected to it in the overall system degrees of freedom. The FEM degrees of freedom and PD degrees of freedom are then coupled through the stiffness contribution of the truss element.

[0102] Furthermore, in this embodiment, the PD near-field neighborhood search and key list construction method in step S103 aims to efficiently construct a flat key data structure that is conducive to GPU parallel scheduling. The specific steps of this method are as follows:

[0103] S131. Accelerated Spatial Grid Construction: To accelerate the neighborhood search process, a uniform spatial background grid is first established in the PD sub-region. All PD material points are assigned to their respective background grid cells according to their coordinates, and a mapping relationship between the material point index and the grid cell is established.

[0104] In this embodiment, a spatial acceleration grid (Cell Grid) is used as an example to illustrate neighborhood search; in other embodiments, other spatial data structures can be used during the S10 host CPU preprocessing stage to accelerate the search of a matter point's near-field neighborhood. All other material points within the grid. These alternative data structures include kd-trees or octrees, which may offer higher search efficiency than uniform grids when dealing with PD models where the distribution of material points is extremely non-uniform.

[0105] S132. Efficient Parallel Search: When traversing each material point i, it is no longer necessary to compare distances with all other material points. Instead, it first locates the background grid cell where material point i is located, and then searches only that cell and all its neighboring material points, calculating the distance between material point i and these candidate material points j. If the distance is less than or equal to the near-field neighborhood radius, the search is successful. If (i, j) is considered to constitute a "key", then (i, j) is considered to constitute a "key".

[0106] S133. Optimized storage of the key list: Store all found "keys" (i, j) in one or more one-dimensional arrays. For example, one array stores the index i of the first material point of the key, and another array stores the index j of the second material point. This "structure array" storage method makes it easy to start N threads in GPU parallel computing steps such as S20 and S40. Each thread independently reads the information of the kth key and performs calculations on it, achieving the best parallel computing efficiency.

[0107] S20. Assembly of the overall stiffness matrix based on CPU and GPU collaboration: This stage aims to construct the overall stiffness matrix K of the system in the initial state. The stiffness contributions of the FEM unit and the interface coupling unit with smaller computational load can be calculated on the CPU, while the stiffness contributions of all "bonds" in the PD region with huge computational load are calculated and assembled in large-scale parallel by starting the GPU kernel function, and finally forming a complete stiffness matrix in the GPU device memory.

[0108] Specifically, in this embodiment, step S20 is a CPU and GPU collaborative overall stiffness matrix assembly method, which aims to efficiently construct the system's overall stiffness matrix K in the initial (damage-free) state. This matrix K is stored in the GPU device's memory and uses a format suitable for parallel sparse matrix operations. This collaborative method allocates computational tasks to the most suitable processor to maximize efficiency, and specifically includes the following steps:

[0109] S201, CPU-side stiffness contribution calculation: On the host (CPU), calculate the stiffness contribution of the FEM sub-region and the interface region.

[0110] FEM subregion stiffness: The CPU iterates through all FEM elements and calculates the element stiffness matrix K for each element. fem .

[0111] Interface coupling stiffness: The CPU iterates through all interface truss elements established in S102 and calculates the global stiffness matrix K of each coupled element according to the method in S123. t .

[0112] Calculate all K fem and K t The index position information of the matrix and its position in the overall stiffness matrix K is stored in a temporary buffer in the CPU memory, ready for subsequent transmission to the GPU.

[0113] S202, GPU-side PD stiffness parallel assembly: On the GPU device, the stiffness contribution of the PD sub-region is calculated and assembled in parallel. This step is achieved by starting the GPU kernel function and allocating a GPU thread for each "key" in the "key list" in S103.

[0114] In this embodiment, a parallel strategy of "one GPU thread corresponding to one PD key" is adopted; in other embodiments, a strategy of "one GPU thread block corresponding to one PD material point" can be used. Under this strategy, all threads within a thread block (e.g., 32 or 64) collaboratively process the key relationships (including stiffness calculation or damage assessment) between the material point and all its neighboring material points. This strategy can better utilize the GPU's shared memory, potentially improving the locality of data access and computational efficiency.

[0115] S203, Stiffness Matrix Integration: Integrating the stiffness contribution data K calculated by the CPU in S201. fem and K t Transfer from CPU host memory to GPU device memory.

[0116] A GPU kernel function is started to atomically add the stiffness contributions from the CPU to the overall stiffness matrix K in the GPU device memory, thus completing the final assembly.

[0117] Furthermore, in this embodiment, the GPU-side PD stiffness parallel assembly method in step S202 parallelizes the PD bond stiffness calculation task, specifically including the following steps:

[0118] S221. Kernel function startup and thread allocation: The host CPU determines the kernel function startup and thread allocation based on the total length N of the "key list". bonds Start a configuration with N bonds The GPU kernel function (Kernel 1) is for each thread. Each GPU thread is uniquely mapped to a PD key. For example, the thread with thread ID k is responsible for processing the k-th key in the "key list" (connecting matter point i and matter point j);

[0119] S222, Parallel Stiffness Calculation: On the GPU, all N bonds Multiple threads execute simultaneously:

[0120] Data reading: Each thread k reads the material point index and coordinates (x, y) of its corresponding key (i, j) from the GPU global memory. i , y i , z i ) and (x j , y j , z j ), and axial stiffness Kp;

[0121] State check: Thread k accesses the "global key state array". In the S20 (first assembly) phase, the state of all keys is "intact" (value is 1).

[0122] Stiffness calculation: If the bond state is "intact", the thread performs stiffness calculation. The interaction between PD material point pairs is treated as a one-dimensional truss element (bond), and its calculation is consistent with the truss element stiffness calculation in S123.

[0123] Calculate the direction cosines l, m, n of the key based on the axial stiffness K. p Using direction cosines, construct the global stiffness matrix k of the bond (i, j). bond .

[0124] S223, Parallel Atomic Assembly: Thread k calculates k... bond The 36 component values ​​are accumulated to the corresponding positions in the overall stiffness matrix K. Since multiple threads may simultaneously write to the same memory address of matrix K (e.g., sharing multiple keys of the same material point i), to avoid "race conditions," this embodiment uses atomic operations. When thread k needs to write to a specific address in matrix K (e.g., the u of material point i), the atomic operation is performed. i When incrementing a value at the diagonal position corresponding to the degree of freedom, it calls the atomicAdd() function. The atomicAdd() operation ensures the atomicity of the "read-modify-write back" operation on that memory address, guaranteeing the accuracy of the concurrent accumulation result. All threads increment their k values ​​in parallel. bond The components of the matrix are safely assembled into the overall stiffness matrix K in the GPU device memory using atomicAdd().

[0125] In this embodiment, the atomicAdd() atomic operation is used to resolve the "write conflict" problem during parallel assembly of the stiffness matrix K. In other embodiments, a graph coloring algorithm can also be used, that is, in the S10 preprocessing stage, all PD bonds are colored to ensure that any two bonds sharing the same material point (i.e., bonds that would cause write conflicts) are assigned different "colors". During assembly in S202, the GPU can execute kernel functions in batches and in color order. Since the memory addresses written by all threads within the same "color" batch do not conflict with each other, the use of atomicAdd() can be completely avoided, and higher assembly performance may be achieved on some GPU architectures.

[0126] S30. Parallel solution of system linear equations based on GPU: This stage is executed entirely on the GPU. Parallelized iterative solvers (such as the conjugate gradient method) are used to quickly solve the large sparse linear equation system Ku=F stored in the GPU device memory, thereby obtaining the displacement vector u of all FE nodes and PD material points in the current system state.

[0127] Specifically, in this embodiment, step S30 is a GPU-based parallel solution method for system linear equations. Its purpose is to quickly solve the large sparse linear equation system Ku=F constructed and stored in the GPU device memory in step S20. This method is executed entirely on the GPU to obtain the global displacement vector u of the current system state. A parallel iterative solver is used, and all core computational steps are mapped to the parallel architecture of the GPU for execution. The specific process includes:

[0128] S301. Initialization: Initialize the vectors required for solving the problem on the GPU, including the initial displacement vector u0 and the residual vector r0 = F - K. u0 And the initial search direction vector p0 = r0.

[0129] S302, Parallel Iterative Solution Loop: An iterative loop is executed on the GPU until the convergence condition is met. Each iteration loop includes the following GPU parallel computations:

[0130] a. Sparse matrix-vector multiplication: parallel computation This operation utilizes a sparse matrix K and a vector p stored in GPU memory. k Multiplication and accumulation operations are performed concurrently by a large number of GPU cores.

[0131] b. Parallel Reduction: The inner product is efficiently computed using a GPU-based parallel reduction algorithm to determine the step size. .

[0132] c. Vector Update: Update the displacement solution vector and residual vector in parallel. Each thread is responsible for one or more components. and .

[0133] d. Convergence criterion: Parallel computation of the norm of the new residual (e.g., through another parallel reduction computation). The iteration is then compared with a preset convergence tolerance. If the result is less than the tolerance, the iteration stops.

[0134] e. Direction Vector Update: If convergence is not achieved, the search direction for the next iteration is calculated and updated through parallel reduction and AXPY operations. (in ).

[0135] S303. Output result: After the iteration converges, the final global displacement vector u is stored in the GPU device memory and can be directly accessed by step S40 without any CPU-GPU data transfer.

[0136] In this embodiment, the conjugate gradient method (CG) is used to achieve parallel iterative solution on the GPU. In other embodiments, other iterative solvers that are also suitable for large-scale parallelization on GPUs can be used to solve the large sparse linear system of equations Ku=F stored in the GPU device memory. These alternative solvers include, but are not limited to, the generalized minimum residual method (GMRES), the biconjugate gradient stabilization method (BiCGSTAB), or the preconditioned conjugate gradient method (PCG). In the PCG scheme, the construction and application steps of the preconditioner (such as the Jacobi preconditioner) are also executed in parallel on the GPU.

[0137] S40. GPU-based parallel computation of damage evolution: This stage aims to determine whether new damage has occurred in the PD region based on the displacement field u obtained in the previous step. By starting the GPU kernel function, the elongation of each "bond" is calculated in parallel and compared with the fracture criterion of the material. If the fracture condition is met, the bond is marked as "fractured" in a global "bond state array" in the GPU device memory.

[0138] Specifically, in this embodiment, step S40 is a GPU-based parallel computation method for damage evolution, which aims to update the material damage state within the PD sub-region in parallel based on the global displacement vector u obtained in S30. This method is implemented by starting a GPU kernel function, and its specific process includes:

[0139] S401, Kernel Function Startup and Thread Allocation: The host CPU determines the kernel function startup and thread allocation based on the total length N of the "key list" in S103. bonds Start a configuration with N bonds The GPU kernel function has multiple threads, and each GPU thread is uniquely mapped to a PD key. For example, the thread with thread ID k is responsible for processing the k-th key (connecting material points i and j).

[0140] S402, Parallel Computation of Bond Elongation: On the GPU, all N bonds Several threads execute concurrently. Each thread k first reads the displacement u of the corresponding material points i and j at both ends of the bond from the global displacement vector u in the GPU device memory. i and u j Subsequently, the thread determines the initial relative positions of the matter points. and current relative displacement Calculate the relative elongation of the bond. .

[0141] S403, Parallel Damage Assessment and State Update: Thread k compares its calculated elongation s with the material's critical elongation s0 (read from GPU memory). If s > s0 and the current state of the key in the "Global Key State Array" is still "intact" (e.g., value 1), then the thread updates the key's state to "fractured" (e.g., value 0).

[0142] S50, Damage Iteration Control and Parallel Update of Stiffness Matrix: This stage is the core of controlling the damage evolution iteration loop. The host CPU first checks whether new damage has occurred in the previous step. If so, the system stiffness changes. At this time, the GPU kernel function will be launched to reassemble the stiffness matrix K in parallel based on the information in the "key state array" (this process will automatically ignore broken keys). After the update is completed, the process will return to step S30 and solve again using the new stiffness matrix. If no new damage has occurred, the convergence of the residual force of the current iteration will be further judged. If the residual force meets the preset threshold, the system is considered to have reached equilibrium and the iteration ends. Otherwise, the process returns to step S30 to continue iterating.

[0143] Specifically, in this embodiment, step S50 is the damage iteration control and stiffness matrix parallel update method, which is the core control process for realizing implicit damage evolution. This process, through the collaboration of CPU and GPU, determines whether the system has reached stability and updates the system stiffness when necessary. Specifically, it includes the following steps:

[0144] S501 Damage Status Check: After the GPU completes the damage assessment in step S40, the host CPU needs to know whether new damage has occurred. This can be achieved by having the GPU perform a parallel reduction operation at the end of S40 to quickly count the number of keys in the "Global Key State Array" whose state has changed (from "intact" to "broken"), thus obtaining a "new damage count". The host CPU then retrieves this "new damage count" value from the GPU.

[0145] S502, Iterative Logic Judgment: The host CPU judges the "new damage count":

[0146] If "New Damage Count" > 0: This indicates that the system stiffness has changed due to material damage, and the current displacement solution u is no longer valid. The system must update the stiffness matrix and solve it again, so the process proceeds to step S503.

[0147] If "New damage count" == 0: and after further convergence judgment of the residual force in the current iteration, the residual force still meets the preset threshold, it indicates that under the current load, the material damage no longer expands, the system has reached equilibrium, the implicit iteration loop ends, and the process jumps to the final S60 step.

[0148] S503, Parallel Stiffness Matrix Update: This step is performed on the GPU and aims to efficiently reconstruct the overall system stiffness matrix K based on the updated "global key state array" in S40, to reflect the stiffness failure of broken keys.

[0149] S504, Loop Return: After the GPU completes the update of the stiffness matrix K in S503, the program control flow returns to S30: the parallel solution step of the system linear equations based on the GPU. The system will use this updated, stiffer-reduced K matrix, combined with the original load vector F, to resolve Ku=F, thereby obtaining the new displacement field u. new After that, the process will re-enter S40, forming a complete "solve-determine-update" iterative loop.

[0150] Furthermore, in this embodiment, the core of the parallel stiffness matrix update method in step S503 is to implement a conditional parallel stiffness assembly, which is specifically implemented as follows:

[0151] S531, Stiffness Matrix Reset: Before reassembly, it is necessary to first clear all stiffness contributions related to PD bonds in the overall stiffness matrix K stored in the GPU device memory.

[0152] S532, Parallel Reassembly Based on Damaged State: The host CPU starts a GPU kernel function (Kernel 3) with similar functionality to (Kernel 1) in S202, and similarly allocates a GPU thread for each key in the "key list". In this kernel function, each thread k (corresponding to key (i,j)) executes a conditional reassembly logic:

[0153] a. Reading the state: Thread k first reads its state value in the "global key state array" from the GPU global memory.

[0154] b. Conditional judgment:

[0155] If the state is "intact" (value 1): thread k performs the exact same operation as S222 and S223, i.e., calculates the global stiffness matrix k of the key. bond Then, the atomicAdd() function is used to safely and atomically add its 36 components to the corresponding positions in the stiffness matrix K.

[0156] If the state is "fractured" (value is 0): thread k will skip all subsequent stiffness calculations and atomicAdd() assembly operations and end the task directly. After this process is completed, the K matrix in the GPU device memory will be updated to a new stiffness matrix that contains only all "intact" keys and contributions from FEM units and interface units.

[0157] In this embodiment, the atomicAdd() atomic operation is used to resolve the "write conflict" problem during parallel assembly of the stiffness matrix K. In other embodiments, a graph coloring algorithm can also be used, that is, in the S10 preprocessing stage, all PD bonds are colored to ensure that any two bonds sharing the same material point (i.e., bonds that would cause write conflicts) are assigned different "colors". During assembly in S503, the GPU can execute kernel functions in batches and in color order. Since the memory addresses written by all threads within the same "color" batch do not conflict with each other, the use of atomicAdd() can be completely avoided, and higher assembly performance may be achieved on some GPU architectures.

[0158] S60. Result Feedback and Post-processing: After the damage iteration cycle ends, the final calculation results, including the displacement field and damage distribution (i.e., the "key state array"), are fed back from the GPU device memory to the host memory for subsequent visualization analysis and data archiving.

[0159] Specifically, in this embodiment, step S60 is the result feedback and post-processing method. It is the final step executed after the iterative logic judgment in S50 indicates that the system has reached equilibrium (i.e., the "new damage count" equals 0). This method aims to copy the final calculation results stored on the GPU back to the host for analysis and archiving. Its specific process includes:

[0160] S601, Data Backhaul: The host CPU initiates a data transfer operation to copy the key result data stored in the GPU device memory back to the CPU host memory. This data mainly includes: the final global displacement vector u (containing the displacements of all FE nodes and PD material points) and the final "global bond state array" (reflecting the "intact" or "broken" state of all bonds in the PD sub-region).

[0161] S602. Post-processing and Output: After receiving the data, the host CPU performs post-processing operations. This includes saving the displacement and damage data to disk files, and calling a visualization library to graphically display the displacement field (generating deformation cloud maps) and the "global key state array" (generating damage cloud maps or crack path maps). Key data points can be extracted for quantitative analysis as needed.

[0162] S603, End: After post-processing is completed, the entire simulation process ends.

[0163] Example 1

[0164] The embodiments of the present invention adopt Figure 6 The failure analysis model of the three-point bending beam with initial cracks is shown below. The geometric dimensions of the component are as follows: Figure 6As shown, the total length is 1000mm, the height D=200mm, the support spacing (span) L=800mm, and the distance from the left and right end supports to the endpoints is 100mm each. A precast vertical crack is set at the bottom of the mid-span of the beam, with an initial crack length a=80mm. The elastic modulus of the material is 30GPa, and Poisson's ratio is 0.3. The fracture behavior of the material is controlled by the critical bond elongation s0 corresponding to its tensile strength of 1.65MPa.

[0165] To apply the coupling method of this invention, the beam is divided into a central near-field dynamics (PD) subdomain and two finite element (FEM) subdomains on either side. The FEM subdomains use 10mm x 10mm four-node rectangular elements, and the particle spacing dx = 2.5mm in the PD subdomains is [missing information - likely a specific value]. =10mm. When solving this model using the S10 to S60 process of this embodiment of the invention, the above-mentioned partitioning discretization, data initialization and transmission to GPU are first completed in step S10. In step S0, the GPU assembles the PD region in parallel (based on the stiffness of dx=2.5mm, the CPU calculates the stiffness of the 10mm×10mm FEM region in collaboration. In the implicit iterative loop of S30 to S50, displacement load is applied. The GPU executes S30 (solving), S40 (damage judgment) and S503 (stiffness update) in parallel to simulate the process of crack initiation and propagation from a=80mm. The GPU-parallel PD-FEM coupled implicit solution method provided by this embodiment of the invention includes the entire simulation process of complex crack propagation. Its core calculations are all completed in parallel on the GPU, which greatly improves the calculation efficiency and can quickly obtain key results such as the load-crack opening displacement (CMOD) curve, verifying the accuracy and efficiency of this embodiment of the invention in processing cracked structures.

[0166] The GPU-parallel PD-FEM coupled implicit solution method described in the above embodiments addresses the shortcomings of existing numerical analysis methods, including the finite element method (FEM) and peridynamics (PD) methods, which suffer from low computational accuracy and efficiency. This invention offers a simple and convenient process, aiming to accurately and efficiently simulate the entire process of solid materials under external loads, from elastic deformation to crack initiation, propagation, and eventual failure. It applies peridynamics (PD) theory to the predicted potential fracture region and the finite element method (FEM) to the remaining linear elastic region, and designs a stable and efficient interface coupling mechanism. Its core innovation lies in addressing the high computational overhead caused by damage evolution in the implicit solution format of the coupled model by proposing a complete GPU parallel computing scheme. This migrates computationally intensive tasks from the CPU to the GPU, fundamentally solving the problem of balancing accuracy and efficiency in fracture simulation in existing technologies, and significantly improving the computational accuracy and efficiency of numerical analysis methods compared to existing technologies. This invention not only achieves a breakthrough in computing speed, but also demonstrates significant advantages in simulation accuracy, numerical stability, and application scope. It provides an efficient, reliable, and powerful numerical simulation solution for solving the long-standing problem of fracture analysis of complex structures in the engineering field.

[0167] Obviously, the embodiments described above are merely preferred embodiments of the present invention, and not all embodiments. The accompanying drawings illustrate preferred embodiments of the present invention, but do not limit the scope of the patent. The present invention can be implemented in many different forms; rather, these embodiments are provided to provide a more thorough and complete understanding of the disclosure of the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing specific embodiments, or make equivalent substitutions for some of the technical features. Any equivalent structures made using the content of this specification and drawings, directly or indirectly applied to other related technical fields, are similarly within the scope of patent protection of this invention.

Claims

1. A GPU-parallel PD-FEM coupled implicit solution method, characterized in that, Includes the following steps: S10. Define the geometric and material properties of the model, discretize it and divide it into PD and FEM sub-regions, construct the bond connection relationship between PD material points, and transfer all the basic data structures required for calculation from the host memory to the GPU device memory. S20. Construct the overall system stiffness matrix K in the initial state. The stiffness contributions of the FEM unit and the interface coupling unit with smaller computational load are calculated on the CPU, while the stiffness contributions of all bonds in the PD region with huge computational load are calculated and assembled in large-scale parallel by starting the GPU kernel function, and finally forming a complete stiffness matrix in the GPU device memory. S30. Use a parallelized iterative solver to solve the large sparse linear equation system Ku=F stored in the GPU device memory to obtain the displacement vector u of all FE nodes and PD material points in the current system state. S40. By starting the GPU kernel function, the elongation of each bond is calculated in parallel and compared with the material's fracture criterion. If the fracture condition is met, the bond is marked as fractured in the global bond state array in the GPU device's memory. S50: The host CPU checks whether any new damage has occurred in S40. If so, the system stiffness changes. At this time, the GPU kernel function is started, and the stiffness matrix K is reassembled in parallel according to the information in the key state array. After the update is completed, the process returns to step S30 and solves the problem again using the new stiffness matrix. If no new damage has occurred, the convergence of the residual force of the current iteration is further judged. If the residual force meets the preset threshold, the system is considered to have reached equilibrium and the iteration ends. Otherwise, the process returns to step S30 to continue iterating and solving the problem. S60. After the damage iteration cycle ends, the final calculation results, including the displacement field and damage distribution, are transferred from the GPU device memory back to the host memory for subsequent visualization analysis and data archiving.

2. The GPU-parallel PD-FEM coupled implicit solution method according to claim 1, characterized in that, The specific steps of S10 are as follows: S101. Define the geometric model, material properties, boundary constraints, and external load conditions of the analysis object. Discretize the geometric model spatially, dividing the model into FEM sub-regions and PD sub-regions. In the FEM sub-region, generate a standard finite element mesh containing nodes and elements. In the PD sub-region, generate discrete material points. S102. Define an interface region at the junction of the FEM subregion and the PD subregion. Within this region, establish the coupling connection between the FEM node and the PD material point. S103. Traverse each material point i in the PD sub-region and search for all other material points in its near-field neighborhood. Store all the searched material point pairs as a global key list. S104. In the host memory, based on the results of discretization and neighborhood search, construct all necessary data arrays, including the coordinate arrays of nodes and material points, the FEM element connection relationship array, the interface coupling element connection relationship array, and the PD bond list generated in S103. In the GPU device memory, allocate contiguous memory space for parallel computing to store the overall stiffness matrix K, displacement vector u, load vector F, global bond state array, material parameters, and coordinates. S105. Perform a one-time data copy operation from host memory to GPU device memory, transferring the coordinate array, all connection relationship arrays, material parameters, and initialized load vector F constructed in S104 from the CPU to the corresponding memory space of the GPU.

3. The GPU-parallel PD-FEM coupled implicit solution method according to claim 2, characterized in that, S102 specifically includes: S121. In the discretized model, explicitly define the FEM sub-region and the PD sub-region, identify the interface located at the junction of the two sub-regions, and determine the FE nodes and PD material points on the interface. S122. Introduce truss elements to bridge the approach field dynamics sub-region and finite element sub-region; S123. Calculate the element stiffness matrix of each truss element using the standard finite element method. .

4. The GPU-parallel PD-FEM coupled implicit solution method according to claim 3, characterized in that, The specific steps of S20 are as follows: S201. On the CPU, calculate the stiffness contribution of the FEM sub-region and the interface region: FEM subregion stiffness: The CPU iterates through all FEM elements and calculates the element stiffness matrix K for each element. fem ; Interface coupling stiffness: The CPU traverses all interface truss elements established in S122 and calculates the global stiffness matrix K of each coupling element according to the method in S123. t ; Calculate all K fem and K t The index position information of the matrix and its position in the overall stiffness matrix K is stored in a temporary buffer in the CPU memory, ready for subsequent transmission to the GPU; S202. On the GPU device, the stiffness contribution of the PD sub-region is calculated and assembled in parallel, which is achieved by starting the GPU kernel function, and a GPU thread is allocated for each key in the key list in S103. S203, The stiffness contribution data K calculated by the CPU in S201 is... fem and K t The stiffness contribution from the CPU is transferred to the GPU device memory, and a GPU kernel function is started. The stiffness contribution from the CPU is then added to the overall stiffness matrix K in the GPU device memory through atomic operations, thus completing the final assembly.

5. The GPU-parallel PD-FEM coupled implicit solution method according to claim 4, characterized in that, The specific steps of S30 are as follows: S301. Initialize the vectors required for solving the problem on the GPU, including the initial displacement vector u0 and the residual vector r0 = F - K. u0 And the initial search direction vector p0 = r0; S302. Execute an iterative loop on the GPU until the convergence condition is met; S303. After the iteration converges, the final global displacement vector u is stored in the GPU device memory.

6. The GPU-parallel PD-FEM coupled implicit solution method according to claim 5, characterized in that, In S302, each iteration loop includes the following GPU parallel computation: Sparse matrix-vector multiplication: Parallel computation This operation utilizes sparse matrices and vectors p stored in GPU memory. k Multiplication and accumulation operations are performed concurrently by a large number of GPU cores; Parallel reduction: The inner product is calculated using a GPU-based parallel reduction algorithm to determine the step size. ; Vector update: The displacement solution vector and residual vector are updated in parallel, with each thread responsible for one or more components. and ; Convergence criterion: The norm of the new residual is calculated in parallel, and then calculated through another parallel reduction. The iteration is then compared with a preset convergence tolerance; if the result is less than the tolerance, the iteration stops. Direction vector update: If convergence is not achieved, the search direction for the next iteration is calculated and updated through parallel reduction and AXPY operations. ,in .

7. The GPU-parallel PD-FEM coupled implicit solution method according to claim 6, characterized in that, The specific steps of S40 are as follows: S401, The host CPU determines the total length N of the key list in S103. bonds Start a configuration with N bonds Each GPU thread is uniquely mapped to a PD key; S402, On the GPU, all N bonds Several threads execute concurrently. Each thread k first reads the displacement u of the corresponding material points i and j at both ends of the bond from the global displacement vector u in the GPU device memory. i and u j Subsequently, the thread determines the initial relative position of the material points. and current relative displacement Calculate the relative elongation of the bond. ; S403. Thread k compares its calculated elongation s with the material's critical elongation s0. If s > s0 and the bond is still in an intact state in the global bond state array, then the thread updates the bond's state to broken.

8. The GPU-parallel PD-FEM coupled implicit solution method according to claim 7, characterized in that, The specific steps of S50 are as follows: S501. The GPU performs a parallel reduction operation to quickly count the number of keys whose states have changed in the global key state array, and obtains a new damage count value. The host CPU retrieves the new damage count value from the GPU. S502, The host CPU judges the newly added damage count: If the new damage count > 0, it indicates that the stiffness of the system has changed due to material damage, and the current displacement solution u is no longer valid. The system must update the stiffness matrix and solve it again, and execute step S503. If the newly added damage count is 0, and after further convergence judgment on the residual force of the current iteration, the residual force still meets the preset threshold, it indicates that under the current load, the material damage no longer expands, the system has reached equilibrium, the implicit iteration loop ends, and step S60 is executed. S503, Parallel update of stiffness matrix to reconstruct the overall system stiffness matrix K; S504. After the GPU completes the update of the stiffness matrix K in S503, the program control flow returns to S30. Based on the GPU-based parallel solution step for the system's linear equations, the system uses the updated K matrix, combined with the original load vector F, to re-solve Ku=F, thereby obtaining the new displacement field u. new And return to S40.

9. The GPU-parallel PD-FEM coupled implicit solution method according to claim 8, wherein step S503 specifically includes: S531. Before reassembly, clear all stiffness contributions related to PD bonds in the overall stiffness matrix K stored in the GPU device memory. S532. The host CPU starts a GPU kernel function and allocates a GPU thread for each key in the key list. In this kernel function, each thread k corresponds to key (i, j) and executes a conditional assembly logic: Reading the state: Thread k first reads its state value from the global key state array in the GPU global memory; Conditional judgment: If the state is intact, thread k calculates the global stiffness matrix k of the key. bond And use the atomicAdd() function to add its 36 components to the corresponding positions of the stiffness matrix K; If the state is fractured, thread k will skip all subsequent stiffness calculations and atomicAdd() assembly operations and terminate the task directly. After this process is completed, the K matrix in the GPU device memory will be updated to a new stiffness matrix that contains only all intact bonds and contributions from FEM units and interface units.

10. The GPU-parallel PD-FEM coupled implicit solution method according to claim 1, wherein the specific steps of S60 are as follows: S601, The host CPU initiates a data transfer operation to copy the key result data stored in the GPU device memory back to the CPU host memory. The data includes the final global displacement vector u and the final global key state array. S602: After receiving the data, the host CPU performs post-processing operations, including saving the displacement and damage data to disk files, and calling the visualization program library to graphically display the displacement field and global key state array. It can also extract key data points for quantitative analysis as needed. S603. After post-processing is completed, the entire simulation process ends.

Images